Ultraviolet cascade in the thermalization of the classical 4 theory in 3 + 1 dimensions. 



in 
o 
o 

(N 
-t— > 

o 
O 

oo 



(N 
> 
O 
OO 
(N 
O 

o 



X 



C. Destri( Q fl and H. J. de Vega( b f[| 

Dipartimento di Fisica G. Occhialini, Universita Milano-Bicocca Piazza della Scienza 3, 
20126 Milano and INFN, sezione di Milano, via Celoria 16, Milano Italia 
(b) lptjje^ Universite Pierre et Marie Curie, Paris VI et Denis Diderot, 
Pans VII, Laboratoire Associe au CNRS UMR 7589, Tour 24, 
5eme. etage, 4> Place Jussieu, 75252 Paris, Cedex 05, France and 
Observatoire de Paris, LERMA. Laboratoire Associe au CNRS UMR 8112. 
61, Avenue de V Observatoire, 75014 Paris, France. 
(Dated: February 2, 2008) 

We investigate the dynamics of thermalization and the approach to equilibrium in the classical 
(f> 4 theory in 3 + 1 spacetime dimensions. The non-equilibrium dynamics is studied by numerically 
solving the equations of motion in a light-cone-like discretization of the model for a broad range of 
initial conditions and energy densities. A smooth cascade of energy towards the ultraviolet is found 
to be the basic mechanism of thermalization. After an initial transient stage, at a time scale of 
several hundreds inverse masses, the squared magnitude of the field spatial gradient becomes larger 
than the nonlinear term and there emerges a stage of universal cascade, independent on the details 
of the initial conditions. As the cascade progresses, the modes with higher wavenumbers, but well 
behind the forefront of the cascade, exhibit weaker and weaker nonlinearities well described by the 
Hartree approximation, while the infrared modes retain strong selfinteractions. As a consequence, 
two timescales for equilibration appear as characteristic of two time-dependent wavenumber regions. 
For k 2 > 02 (t), we observe an effective equilibration to a time-dependent power-like spectrum 
with a time scale in the hundreds of inverse masses; cutoff effects are absent and the Hartree 
approximation holds for k 2 3> 4> 2 (t). On the other hand, infrared modes with k 2 < 4> 2 (t) equilibrate 
only by time scales in the millions of inverse masses when cutoff effect become dominant and 
complete thermalization is setting in. Accordingly, we observe in the field correlator a relatively large 
and long-lived wavefunction renormalization of nonperturbative character. There corresponds an 
effective mass governing the long distance behaviour of the correlator which turns to be significantly 
smaller than the Hartree mass which is exhibited by the modes with k 2 > (f> 2 (t). Virialization and 
the equation of state start to set much earlier than thermalization. The applicability of these results 
in quantum field theory for large occupation numbers and small coupling is analyzed. 



Contents 



I I. Introduction! 

I II. Th e model in the continuum a nd on the lattice! 



B. Rcfiularization on a finite s pace lattice ] 
C^^iscretizec^d^iamics on a spacetime lattice! 
j^^Weragecl^bserrable^ 



IlII. Th ermal equilibr ium! 
A^ 



B. Perturbation theory. Hartree rcsummation and beyond I 



C. ETm^ibrk^i^ri^^clatTicel 
D"^oTrt^^arl^^miul^tronsl 



llV. Dy namics of Th ermalization! 

A^mTraj^OTurrtionsI 



B. Time evolution of local observablcs 



C. Virialization and Equation of State 



4 
4 
■5 
G 

8 

9 
9 

11 
12 
14 

16 
16 
17 

20 



'Electronic address: Claudio.Destri@mib.infn.it 
^Electronic address: devega@lpthe.jussieu.fr 



2 



D^Jj^yahrti^^ 23 
^ 29 



^ 32 



G. Effective frequency and particle numbed 33 

H. Effective Mass sauared| 38 



I. The slow dynamics of the infrared modes and effective thermalization! 39 



V. Discussions and Conclusions! 40 

A. Average over discrete directions! 42 

B. Linearized discrete dynamics! 43 

C. The Hartree approximation from the late time exact behaviour! 45 

D. The cnoidal solution! 46 
References! 46 



I. INTRODUCTION 



The understanding of the dynamics of thermalization and relaxation in a field theory is a subject of critical impor- 
tance both in early cosmology as well as in ultrarelativistic heavy ion collisions. 

Pioneering work in this topic was initiated by Fermi, Pasta and Ulam[lJ for a chain of coupled oscillators. Since 
then, this problem has been studied within a variety of models 2] with the goal of answering fundamental questions 
on ergodicity, equipartition and in general the approach to equilibrium in non-linear theories with a large but finite 
number of degrees of freedom. 

By understanding the dynamics of thermalization we mean to uncover in detail the mechanisms that leads to 
thermal states starting from arbitrary initial conditions as well as to determine the time scales of these phenomena. 
In cosmology the inflationary paradigm assumes that after the inflationary stage a period of particle production 
and relaxation leads to a state of local thermal equilibrium thus merging inflation with the standard hot big bang 
cosmology 00- 

Inflationary scenarios lead to particle production either via parametric amplification of fluctuations of the inflaton 
(in the case of an oscillating inflaton) or by spinodal instabilities during phase transitions 0,0, 8,9]. In both cases the 
non-equilibrium dynamics is non-perturbative and results in a large population of soft quanta whose dynamics is nearly 
classical. The non-equilibrium dynamics of particle production and eventual thermalization are non-perturbative and 
the resulting fluctuations contribute to the evolution of the scale factor^amely the back- reaction from the fluctuations 
becomes important in the evolution of the cosmological space-time 0, 0, H @- Both parametric amplification or 
spinodal decomposition lead to non-perturbative particle production of a band of wave-vectors, typically for soft 
momentajl 0, H 0]. This non-perturbatively large population allows a classical treatment of the non-equilibrium 
evolution. Thermalization may be described in this classical framework or one can implement the quantum 2PI 
framework in cosmological spacetimes. 

In this article we study the non-equilibrium dynamics going towards thermalization in the classical </> 4 theory in 
3 + 1 dimensions. As mentioned above the initial stages of non-equilibrium dynamics either in cosmology or in 
ultrarelativistic heavy ion collisions is mainly classical. Classical field theory must be understood with an ultraviolet 
cutoff, or equivalently an underlying lattice spacing, to avoid the Rayleigh- Jeans catastrophe. 

The advantage of studying the classical field theory is that the equation of motion can be solved exactly. Contrary 
to quantum theory, there is no need of making approximations. Classical field theory is expected to be a good 
approximation to the quantum theory for large occupation numbers and small coupling as we argue in this paper. 

We provide a detailed understanding of the main dynamical mechanisms that lead to thermalization and to study 
the approach to thermalization by several different observables. We address several questions on the equilibrium 
and non-equilibrium aspects: i) what are the criteria for thermalization in an interacting theory?, ii) what is the 
mechanism that leads to thermalization?, iii) what is the dynamics for different observables, how they reach thermal 
equilibrium and which are the relevant time scales?. 

We focus our study on these issues within the context of classical field theory, which is an interesting and timely 
problem all by itself. It is also likely to describe the initial stages of evolution strongly out of equilibrium in quantum 
theory for the relevant cases of very large occupation numbers and small coupling. 
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The main results of this work are the following: 

• We implement a discretization of the p 4 model in any number of space-time dimensions which is very accurate 
and stable, maintains the relativistic symmetry between space and time and conserves energy exactly (that is 
to machine accuracy on a computer). The same discretization scheme, restricted to space alone, can be used to 
study the model at canonical equilibrium through Monte Carlo simulations. 

• We provide the main physical properties of the classical ip model in thermal equilibrium using low temperature 
perturbative expansions as well as Monte Carlo simulations. These simulations show that the leading order 
perturbative results have a large domain of validity, further extended by the Hartree approximation, except for 
infrared modes, which retain a strong nonperturbative character. 

• We extensively investigate the dynamics of the ip model for a wide range of (infrared supported) initial conditions 
and energy densities. After a first stage with relatively important fluctuations whose precise structure depends 
on the details of the initial conditions, the model evolves towards a universal stage where energy transfer 
from low to high wavenumbers becomes steady and very effective, resulting in a steady smooth ultraviolet 
cascade. Namely, the power spectrum of the field (f> and its canonical momentum 7r acquire support in a 
monotonic and slow way over larger and larger values of the wavevectors. This ultraviolet cascade leads to a 
efficient transfer of energy from ir 2 , <fi 2 and the interaction terms <p 4 to the spatial gradient (V0) 2 which grows 
monotonically. Therefore, there is a crossover during this stage from a strongly to a weakly interacting theory 
since the non-linear term 4> 4 becomes much smaller than the spatial gradients. The average wavenumber k(t) of 
the modes grows monotonically with time approximately as i 1 / 3 . The universal stage is well established by a 
time to ~ 500 (in inverse mass units) which does not depend on the lattice spacing and depends weakly on the 
energy density E/V as long as E/V is not very small. The study of the <j> and ir = <fi correlators allows us to 
obtain the power spectrum of <j> and n which exhibits universal scaling properties during the smooth cascade. We 
find that virialization sets rather fast at times t > to- Parallel to virialization, the equation of state approaches 
the radiation equation. Namely, the ratio of the pressure divided by the energy density approaches | for energy 
densities not too small as the nonlinearities weaken. The ultraviolet cascade continues till the lattice cutoff is 
reached and therefore the universal properties w.r.t. the discretization method are lost. These cutoff effects 
become important long before full thermalization of the power spectra sets in, but late enough, for small enough 
lattice spacing, for the window of the universal cascade to be clearly visible. 

• The infrared modes with k 2 < 4> 2 {i) exhibit a much slower dynamics than the rest of the modes. The modes with 
kit) 3> k 2 3> (j) 2 are well described by the Hartree approximation and equilibrate effectively by times ~ to — 500 
(in inverse mass units) while the infrared modes reach an equilibrium state only by times t ~ 10 6 , when k(t) is 
very close to the UV cutoff. In this equilibrium state a rather large wavefunction renormalization remains for 
low wavenumbers. In configuration space the correlation length turns to be 1/M c g(t) where the effective mass 
of the infrared modes M e g(t) is substantially smaller than the Hartree mass. Therefore, there are two scales 
for equilibration: the shorter one to characterizes the UV cascade evolution, the local magnitudes as the time 
average <fi 2 (t) and the modes with k 2 > </> 2 (t) , while the longer scale governs the evolution of the infrared modes 
and the wavefunction renormalization for k 2 < 4> 2 (t). 

• The thermalization process in 3 + 1 dimensions is quite different from that in 1 + 1, studied with the same light- 
cone-like discretization approach in ref.Q, although a universal smooth ultraviolet cascade is present in both 
cases. In the 1 + 1 case the cascade of the n power spectrum is characterized by a single universal shape function 
with the time evolution reducing to a scale transformation on such function. On the contrary, in D = 3 + 1, one 
can still define a shape function for the cascade but it turns out to be time-dependent; in particular, there exist 
a window at intermediate values of the scaled wavenumber k/k(t) with a power-like behaviour, but the power 
depends markedly on time even when cutoff effects are fully negligible. Most remarkably, the shape function is 
almost flat near the origin in 1 + 1 dimensions, allowing to consistently define an effective thermalization which, 
starting from the deep infrared where power is concentrated at zero time, progresses to higher wavenumbers with 
a well defined effective temperature monotonically decreasing in time. In D = 3 + 1 instead, the deep infrared 
is the last to thermalize and the mode equilibration in the bulk of the cascade, where power-like spectra are 
observed, cannot be characterized only by a time-dependent effective temperature. 

Nonetheless, the notion of a time-dependent effective temperature keeps its sense for specific (simple) observables 
also in 3 + 1 dimensions. For instance, from the late time behaviour of ir 2 , cf> 2 , <fi 4 and (V</>) one can read 
off an effective temperature which monotonically decreases with time approximately as ~ i". This effective 
temperature reaches the proper nonzero limit in the lattice for very late times showing that true thermalization 
has been achieved. In the continuum this effective temperature would always vanish for infinite time. 
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• Although all our present work (and in ref. 6j) deals with classical field theory, some of the main features as the 
ultraviolet cascade and the slow thermalization dynamics of the infrared modes should be relevant in quantum 
field theory for large occupation numbers and small coupling. 

The main message to early cosmology and ultrarelativistic heavy ion collision physics from the present work is 
that thermalization proceeds very slowly in classical field theory (at least for unbroken symmetry in scalar models). 
However, a state of effective equilibration is soon reached and is characterized by a steady, smooth and relatively 
simple cascade of energy flowing towards the ultraviolet. The natural interpretation is that the "fast" degrees of 
freedom have indeed reached some sort of statistical equilibrium defined by few macroscopic parameters (such as the 
temperature, but not the temperature alone) which are slowly varying in time. On the continuum, without any cutoff 
effect for arbitrary long time, this effective temperature should become the only relevant slow-varying quantity for 
late enough time; it would eventually vanish for infinite time in classical theory, while in quantum theory it should 
reach a finite nonzero value necessarily involving h. 

Thermalization has been reached in quantum field theories by numerical studies in 3 + 1 and 2 + 1 dimensions in refs. 
[Til ITU and 01 1 respectively. In these works the 2PI expansion to lowest order and the analogous Kadanoff-Baym 
approximation are used, respectively. The couplings considered are quite strong and do not allow a direct connection 
with the classical dynamics investigated in the present paper. The numerical studies previously reported on the 
dynamics of classical and quantum field theories |l9j]-|23j| have not yet focused on studying the mechanism of energy 
transfer from long to short wavelengths as a function of time. In ref. 19] the emergence of spatio-temporal structures 
is reported for a symmetry breaking (j) 4 model. Refs. [2(J study the evolution of the classical 1 + 1-dimensional 4> 4 
model. For further work in this domain see refs. [T^. IT3| and |28l |. 

It would indeed be interesting to study if the universal cascade found in the classical theoryQ remains at least 
during some early and intermediate stages in the quantum theory. 

Ref. ^5 investigates the classical in the 3 + 1-dimensional massless <ft 4 model. The ultraviolet cascade is present in 
their numerical results. In addition, an analytic scenario for the scaling behaviour is proposed following the transport 
treatment of weak wave turbulence |25|. However, this heuristic treatment does not fully capture the thermalization 
dynamics found in 1 + 1 dimensions in |6| as well as in 3 + 1 dimensions in the present paper. 

This paper is organized as follows: in sec. II we present the classical (f> 4 model in the continuum and in the 
lattice and our way to perform suitable coarse-grainings. Sec. Ill discuss the classical <fi 4 model in thermodynamic 
equilibrium and its properties both from perturbative as well as Monte Carlo calculations. Sec. IV is the core of the 
article where the dynamics of thermalization is presented. We study the time evolution of local observables as well as 
the correlators whose Fourier transforms yield the power spectra of the field <fi and its conjugate momentum it. The 
early virialization and the late thermalization of the infrared modes is discussed. Sec. V contains discussions and 
conclusions while four appendices deal with more specific topics. 

II. THE MODEL IN THE CONTINUUM AND ON THE LATTICE 

Besides defining the classical (f> 4 model, we present in this section its lattice version simultaneously discretizing 
space and time in such a way that an exactly conserved lattice energy can be defined. The average procedure over 
the basic physical observables is presented too. 

A. Basic definitions and notations 

The Lagrangian density of the (3 + 1)— dimensional ip 4 field theory reads 

£ m ,A m = -d^ipd^cp - — V - ^ V 
where = d/dx^, \i = 0, 1,2, 3. This leads to the classical equation of motion 

a^ + m 2 ^ + A^=0 (2.1) 

In 3 + 1 space-time dimensions and standard units H = c = 1, the coupling A is dimensionless while <p has the 
dimensions of a mass. 



5 



At the classical level one can always rescale the coordinates and the field using some reference mass M and absorb 
the coupling A in the field. Thus setting 



M 

<p(x) = -j= <f>(Mx) 



(2.2) 



and renaming (Mt, M x) as (t,x), yields for the dimensionless field 4> the equation 

4>- V 2 + 70 + 3 =0 

where 4> = d<j>/dt, V .,■</> = d<j>/dxj and 7 = m 2 /M 2 . If m? > 0, then we can choose M = m and study the 
parameter-free equation 







(2.3) 



This is our choice, having assumed a massive (p. At any rate, the important point here is that, unlike in QFT, the 
notion of coupling constant at the classical level is not absolute as it can be scaled out. 

In terms of <p ancl its canonical conjugate momentum 7r = </>, the Hamiltonian is the standard sum of a kinetic plus 
a potential term 



T[ir] +V[4>], T[tt] 



d 3 xir 2 , V[4>] 



1 



d 3 x 



(V^) 2 + 2 + \ 4 



(2.4) 



This Hamiltonian is dimensionless; the energy of the original field ip in the standard dimension-full coordinates is 
given by (m/X) H[tt, (j>\. 

We consider the model restricted to a finite volume, which we take to be the cube of side L (in units of to -1 ) 
and volume V = L 3 ; we assume periodic boundary conditions (p.b.c), namely <f)(x + Ln,t) = <fi(x,t) for any 
n = (ni, 712,713) G Z 3 . 

In terms of standard Fourier mode amplitudes, 



TTfc 



d 3 



X 



^ — ik-: 



the Hamiltonian reads 



d 3 x 

v 1 / 2 



—ih-; 



(x) 



-k 1 



k = — n 

Li 



n e 1? 



(2.5) 



H l<f>, *] = \ E [l^l 2 + 0- + k 2 Mk\ 2 + 2V E 1 



>g <j>q' <Pk 9-q-q'-k 



k qq' 

The wavenumbers k are dimensionless; their dimension-full counterpart are given by mk. 



B. Regularization on a finite space lattice 

For any numerical treatment it is necessary to discretize space and time over some lattice. This obviously introduce 
an ultraviolet cutoff A over the wavenumbers of the order of the inverse of the lattice spacing. To discretize the space 
we consider the cubic lattice (2a Z) 3 , a being half the lattice spacing. Taking into account the finite size L, we obtain 
regularized fields 

4> n = 4>(x) , n n = n(x) , x = 2an, n e Cn , 

where Cn is the cubic subset of Z 3 formed by integer triples (m, ri2, n^) satisfying —N/2 + 1 < rij < N/2, j = 1, 2, 3. 
Here N = L/ (2a) is assumed even and N 3 is clearly the total number of degrees of freedom of the regularized fields. 

Owing to p.b.c, to the cube 2a CV in 2;— space there corresponds the dual cube (2tt/L)Cn (the so-called first 
Brillouin zone) of wavenumbers k fulfilling 

2tt 

k=—n, -N/2+ 1 < nj < N/2 , j = 1,2,3 . 
L 

The largest value of each component of k, A — (2n/L)(N/2) = n/(2a), is the UV cutoff. 

The lattice form of the Hamiltonian is not unique, being restricted solely by the requirement to formally reduce to 
the continuum expression eq. (|2.4|) in the a — > limit. For all ultralocal terms, with fields at coincident points one 
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could assume the simplest discretization as sums of one-site terms. Thus the kinetic energy T would read on the 
lattice 



2 

and similarly for the space integrals of (f> 2 and <f> . The integral of the gradient term (V0) 2 is converted, through 
integration by parts, to the integral of — </>V 2 </>, with V 2 the Laplacian. Then, in order to define the theory on the 
lattice, V 2 is replaced by a discretized Laplacian. The simplest choice is the nearest-neighbor form that leads to the 
well known replacement of the spectrum fc 2 on the continuum 

2 j. 2 » _ sin ka 
a 

However, our choice of discretization is different, as we show in the next subsection. 



C. Discretized dynamics on a spacetime lattice 



In order to solve numerically the evolution equations for the 4 theory, it is necessary to discretize time besides space. 
In most approaches, space and time discretization arc performed separately. As discussed in the previous section, space 
discretization turns the field theory into a classical dynamical problem with finitely many 'coordinates' and 'momenta'. 
Their evolution is governed by ordinary differential equations which eventually require some discretization of time 
to be solved numerically. Here we proceed in a different way, treating space and time in a symmetric way. This 
generalizes the scheme introduced in 1 + 1 dimensions 0, |2(| . 

In our approach space and time are simultaneously discretized, with the same lattice spacing 2a, in a staggered fash- 
ion over the lattice Z D+1 U (Z + 1 /2) D+1 (for the sake of generality, we consider here a space of generic dimensionality 
D). In other words the discretized space-time points are 

(x, t) = a (n, s) 

where the integer components of (n, s) are either all even or all odd. This allows, as we shall see, to protect to a large 
extent the relativistic invariance of the continuum field equation (|2.3(l and to provide an exactly conserved energy on 
the lattice. 

Let us first of all define the averages over local cubes of field powers 



where cr = (ai, er 2 , . . . , a D ), a l = ±, i = 
observables through the correspondences 

4>(x, t) = tt(x, t) 
4> 2 {x,t) 
^(x,t) 
(V<t>) 2 (x,t) 



a*,t)Y>, p = 0,1,2,... (2.7) 
1,2, ... ,D. We then construct two lattice versions of relevant continuum 



n±(x,t) = ±a~ 1 [4>(x,t±a) - <$> {1 \x,t)] 
4> 2 ± (x,t) = \ [<P 2 (x,t±a) + $ (2) (x, t)] 
</>±{x,t) = <p 2 (x,t±a) $ (2) (x,t) 
(D</>) 2 (;r,i) ee a- 2 {& 2 \x,t) - t)] 2 } 



By construction, the lattice quantities on the r.h.s. tend to the continuum expressions in the limit a 
the lattice energy densities 



(2.8) 



0. Therefore 



1 



£ ± (x,t) = -(2a) D [Tr 2 ± + (D<t>) 



4 



both tend to the continuum energy density as a —> 



£±~£ = 



1 



d D x 
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upon the natural identification (2a) D ~ dPx. Notice that £± can be explicitly written as 

r -| 2 

£±(x, t) = \a D ~ 2 0(a;, i ± a) — (f>(x + cut, t) 



+ i (2a) D [l + </> 2 (x , t ± a)] 1 + -L £ 2 ( x + aCT , t) 



i(2«) D 



(2.9) 



which shows that only diagonal, space-time symmetric finite differences are present. 
To higher orders in a, £ and £± do differ; in fact 



£+(x, t) - £-(x, t)=4 (2a) u - 2 [<j>(x, t + a) - cf>(x, t - a)} Q(x, t) 



where 



Q{x,t) = \ [(/>(x,t + a) + <l>(x,t- a)]|l + \a 2 [l + <f> 2 {x, t)] } - $i(cc,<) 
Hence, if Q(a?, i) = 0, then £+(x, t) = £-(x, t) also on the lattice. In this case the total lattice energy 

E= Yl e +m 



x£2ai 



is exactly conserved in time, since it can also be written 



E= ^ £-{x 1 t + a) 

x£2a(Z+l/2) D 



(2.10) 



(2.11) 



All this holds exactly on infinite space. If space is restricted to the (hyper)cube of side L, suitable boundary conditions 
on </>(x, t) are necessary; p.b.c. are of this type if L = 2Na with N an integer. 

In conclusion, we may regard Q(x, t) = as a discrete field equation which conserves the total energy E. Explicitly, 
Q(x, t) = reads 



(x, t + a) + 4>(x, t — a) = 



2^E CT <t>{x + a<r,t) 



l + \a 2 



l + 2- D E CT <t> 2 {x + acT,t) 



(2.12) 



which evidently allows to evolve in time any configuration known on two consecutive time slices, say t = and t = a. 

It is easy to check that Q(x,t) = indeed becomes eq. (|2.3|l in the continuum a — ► limit. The order a is trivially 
satisfied, odd powers of a vanish identically as a consequence of the symmetry of eq. I|2.12|l under a — > —a, while the 
order a 2 produces eq. (|2.3fl . 

Keeping up to C(a 4 ) in ea. (|2.12|) yields, 



V 2 + , 



a 2 Q 2 + 0(a 4 



2 * / v 12 v t 4 



2v72, 



v z v 



2 ^ ay 



[(V0) 2 + cb V 2 cf\ 



(2.13) 



To cast ea. H2.120 in a form suitable for numerical calculations, we define the lattice arrays F(n, s) and G(n, s) as 

F(n, s) = 4>(2na, 2sa) , G(n, s) = 0(2na - er + a, 2sa + a) , (2-14) 



where n G Z 15 , seZ and cr + = (1,1,,..., 1). 
We then obtain the iterative system 

F(n,s + 1) = -F(n,s) + 

G{n, s + 1) = -G(n,s) + 



2^-i + ±a 2 [2^ + E T G 2 (n + r,s)] 



E T %-r,s + l) 



2D-1 + i a 2 [2£> + £ t F 2 {n -t,s + 1)] 



(2.15) 



where r = (ti, r2, . . . , 7~d), with r, = 0, 1 and, according to the p.b.c, F(n + Nt, s) — F(n, s) and G{n + Nt, s) 
G{n, s) for any r. 
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As initial conditions we have to specify F(n,0) and G(n,0) for < rij < N — 1, i = 1,2, ... ,13. Once these 
values of the fields are provided, the iteration rules eqs. I|2.15|l uniquely define F(n, s) and G(n, s) for all s > 0. 
A comparison of this discretized dynamics with other more traditional numerical treatments of hyperbolic partial 
differential equations was performed in 1 + 1 dimensions Here we only notice that this approach is particularly 
efficient, stable and accurate, specially when the continuum limit a — > and very long evolution times are of interest. 

All observables of the continuum can be rewritten on the lattice in terms of the basic fields F(n, s) and G(n, s), 
according to the correspondences rules eqs. (|2.8fl and the identification eq. (|2.14|) . In the sequel, while referring to 
observables discretized as above, whenever possible we shall keep using the notation corresponding to continuum 
observables for simplicity. Particular care must be taken for </> 4 , since its lattice definition in eqs. I|2.8[> uses a product 
of fields over different lattice sizes. This is harmless in the continuum limit for smooth fields, but makes a difference 
for fields with large ultraviolet support, as we shall see shortly. For this reason, we shall keep the notation (f)±{x,t) 
given in eqs. (|2.8|) well distinct from the ultralocal definition 4 (a;,t) = [4>(x 1 t)] A . 

D. Averaged observables 

The key observables in our investigation are the basic quantities 

0, <f> 2 , 4 , 2 , (V0) 2 , (2.16) 

as well as the power spectra of <f> and tt, that is |</>fc(i)| 2 and |7Tfc(i)| 2 , where <fik(t) and 7Tfc(t) are given by ea. H2.5|) . Let 
us recall again that we are using here the continuum notation also for the Fourier transforms, although they actually 
are discrete Fourier transforms. 

The fluctuations of all these observables do not vanish upon time evolution. Hence for generic initial conditions 
they do not have any limit as t — > oo. These are fine-grained or microscopic observables. Typically there are several 
spatio-temporal scales, the microscopic scales correspond to very fast oscillations and short distance variations that 
are of no relevance to a thermodynamic description. We are interested in longer, macroscopic scales that describe the 
relaxation of observables towards a state of equilibrium. 

In particular, the ergodic postulate states that ensemble averages must be identified with long time averages as 
well as spatial averages over macroscopic-sized regions. Hence to make contact between the time evolution and the 
thermal averages we need to properly average the microscopic fluctuations. 

First of all, for local quantities such as those in ea. (|2.16|) we take the spatial average. Secondly, we take suitable 
time averages of all key observables in the following way 

d?(t)=- [ dt' - [ d 3 ^ 2 (a;,t'). (2.17) 

T Jt-r y JV 

where r 3> a fixes a limit to our resolution power in time, r need not be constant throughout the time evolution. For 
instance, we find that a practical and efficient choice is 

T(t)=T(0) + Ct (2.18) 

where C is small and positive with typical values ~ 0.1. In this way we still keep t r(t) while the dependence on 
the initial values becomes negligible for practically accessible times. This method is quite effective in revealing general 
features of the (logarithmic) time evolution such as the presence of distinct stages characterized by well separated 
macroscopic time scales. 

Besides the time average, for more detailed observables such us the power spectra we performed an average over 
the discrete directions in space as discussed in Appendix^] Moreover, we sometimes averaged also over several initial 
conditions. All together, we denote the results of all these coarse-grainings simply with an overbar to avoid cluttering 
of notation. For example, assuming equiprobable initial conditions 

4^!f/4 i d3x few)]® • 

where the superscript W labels the M different choices of smooth initial conditions, all with the same energy density 
E/V. Likewise, recalling the wave-vectors quantization k = (2n/L) n, we have 

W(t) = TrE ^^^(n<|n|<n + l)- / dt' \^\t')\ 2 (2.19) 

M i= l bn n T Jt-r 
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where the discrete radial wavenumber k is the average |fe| over the shell of all S n wave- vectors k = (2ir/L) n satisfying 
n < \n\ < n + 1 [see Appendix 1X1 for more detail]. Analogous expressions hold for </>, <j) 2 , </> 4 , tt 2 , (V</>) 2 , |7Tfc| 2 and 

\M 2 - 

In particular, due to the linearity of these averages, we have the sum rules: 

+A ,31, _^ r +A j3jL 



(<)=<«') , / — fel 2 (t) = * 2 (t) , (2.20) 



where A = 7r/(2a) is the UV cutoff on the lattice and we write the wavenumbers as continuous although they are 
discrete in all actual calculations. 

Finally, let us recall that power spectra are just the Fourier transforms of space-averaged correlation functions. For 
instance: 



\4>u\ 2 {t) = / fx e-*" <M>(x,t) , (2.21) 
Jv 

where, according to our general rules, 



M 



III. THERMAL EQUILIBRIUM 

We present in this section the basic properties of the <j) A theory in thermal equilibrium. These results will be 
compared in the subsequent sections with the time averages in order to asses whether and how thermalization is 
approached. 



A. General aspects 

The thermal average of any physical quantity O = Q[cj>, tt] in the canonical ensemble is written as 

where / J Dcj) Dir stands for functional integration over the classical phase space and j3 = 1/T is the inverse (dimcn- 
sionless) temperature in the dimensionless variables. In terms of the physical temperature, here defined as T p , T is 
given by 



(3 m 



T=- = -T p . (3.2) 



As a consequence of the field redefinition available in the classical theory, the relevant variable for equilibrium ther- 
modynamics is T. Therefore, for a fixed physical temperature T p we see from eq. (|3.2|) that the low temperature limit 
T <C 1 corresponds to the weak coupling limit and/or T p <C m. This will be relevant in the analysis below. 

Translation invariance (which is preserved by p.b.c.) implies that averages of local observables Q(x) which depend 
on tt and <f> only at one point x, do not depend on x, that is (Q(x)) — (6(0)) = (O). 

Furthermore, the fact that the Hamiltonian is the sum of a kinetic and a potential term [see Eq. (|2.4|l ] entails that 
the average of observables which are of the form 

e[<M = ei[(fl e 2 [7r] , 

factorize as 

(6[0,7r]) = (9i[0]) (6 2 M). 

with 
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Moreover, since the 7r— integration is Gaussian and ultralocal, it can be performed quite easily in most cases, leaving 
the configurational integral over <fr = 4>{ x ) for x <= V to be computed. 

It follows from eqs.JQJ and (|3.1[l that the two-point correlation function of the canonical momentum tt(x) in the 
classical theory in equilibrium is given by 

(n(x) tt(x ')> = T8(x - x') . (3.3) 

which leads to a flat power spectrum for -k 

<l^| 2 >=7\ (3.4) 

where we used the Fourier transform eq. (|2.5|) . This is of course a consequence of equipartition, and gives a criterion 
to identify the temperature: the height of the flat region in the power spectrum of 7r if such flat region shows up. 

Eq. I|3.4[l implies that (it 2 ) = T 6(0) where 6(0) is to be understood as made finite by some UV cutoff procedure. If a 
rotationally invariant sharp cutoff A is used, then £(0) = A 3 / (6n 2 ). In case of the lattice regularization of section lll Bl 
with ir(x) entering only in the kinetic energy converted to a sum as in eq. (|2.6I) . we have instead 6(0) = l/(2a) 3 . 
Thus, 

^(sharpA) = , 2) (lattice. 2a) = T (35) 

6 7r (2a) 

Equating these two regularized averages would yield the the identification A = (6/tt) 1 ^ 3 (2tt/o) = (6/7r) 1 / 3 A; this 
connection between cutoffs is however not universal; it is valid only for (ir 2 ) and in general different for other UV- 
divergent quantities. Moreover, we will see below that the discretized dynamics described in section HTCl provides a 
slightly different UV regularization even for (|7Tfc| 2 ). 

At equilibrium, the classical virial theorem takes the form 

(0 2 ) = ((V0) 2 ) + (cf) + (0 4 ) . (3.6) 

where we have trivially generalized to three spatial dimensions the derivation in ref. [6j. When combined with the 
energy functional H [-7T, <p] given in eq. it yields 



v =v = {7r) ~* 



^>-i(0 4 >, (3.7) 



where E is the average energy in the canonical ensemble. Since N = V/ (2a) 3 is the total number of degrees of freedom 
in the simplest lattice regularization, we find using ea. H3.5f) that the temperature T is related to the (average) energy 
per degree of freedom as 

T=| + 2a 3 (0 4 ). (3.8) 

Therefore, for a <C 1, close to the continuum limit where one finds a 2 (4> 4 ) to be finite as a — > (see below), the 
temperature is identified with the energy per site. Furthermore, 

T=(2a) 3 [p+i(0 4 )] , (3.9) 

where p = E/V is the energy density, while the pressure p is given in general (not necessarily at thermal equilibrium) 
by 



2_ 1 ( V0) 2 _02_ 1 04j _ (310) 

At thermal equilibrium we find from eas. H3.6fl . (|3.9|l and (|3.10|) the equation of state 

(p) = L p -^ 2 ) (3.11) 

that is a radiation-dominated equation of state (p) ~ p/3, since ((f) 2 ) is of order o , much smaller than (p) and p 
which arc both of order a -3 for fixed T. 
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B. Perturbation theory, Hartree resummation and beyond 

The continuum configurational functional integrals 

can be computed in a perturbative expansion in powers of T. In order to do that one changes the functional integration 
variable 4>{x) to 

X{x) = -1= <t>{x) . 
v 2 



Ea. (|3.12|l takes thus the form 

<*'> = ' j^—n^w^ <3 - 13 > 

We read from ea. (|3.13[l the associate Feynman rules. The Euclidean \ propagator in three space dimensions takes in 
momentum space the form 

A 1 



k 2 + 1 



and the quadrilinear \ vertex has (—6 T) as coefficient. 

As already recalled above, this defines a superrenormalizable field theory with only two divergent diagrams: the 
one-loop tadpole and the two-loops sunset. Recall that in classical statistical mechanics the regularized (bare) theory 
is the physical one. The ultraviolet divergences are physical and cannot be eliminated by counter-terms as in quantum 
field theory. 

The two-point function ((f>(x)(j)(x / )) can be written as 

/rfik i T 

^<|fc|'>e*<") , (l^l 2 )= fc2 + 1 + s(fc2) (3-14) 

where S(fc 2 ) stands for the self-energy. 

To lowest order in T the tadpole takes then the form 

W^/.^^^IHOIA-)] (3.15, 



t<s (2wf B + 1 2k 



while the self-energy reads to lowest order 



S = 3(^) A ^ 1 ^[1 + 0(A- 1 )] 

Here we are using the sharp spherically symmetric UV cutoff, to be compared later on to the lattice regularization 
method with cutoff A = it/ (2a). 

To next order £ contains the two- loops sunset diagram, which is only logarithmically divergent; therefore the 
tadpole dominates the effective mass in the limit A — > oo. Hence, in this limit the two-point function is dominated 
by the Hartree approximation (the sum of all daisy diagrams) 



(I ^ |2 ^FTTTW) (3 ' 16) 



where the self-consistent (d> 2 ) fulfills 



d 3 k „- f d 3 k T T r 



(l^r) = 7^17, I i , o/^v A-\ ITA+... (3.17) 



(2tt) 3 xiy ^' ' J (2tt) 3 k 2 + 1 + 3 ((/) 2 ) ~ 2ir 2 
We analogously find for the gradient squared (V</>) 2 , at lowest order 



<(V0) 2 ) = ^ (A 2 - Aa ) +C)(.V'!. ( 3 .IS) 
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and for (<f) A ) using Wick's theorem, 



6 2 ) 2 [1 + OiA- 1 )] = -^j- [A 2 + 0(A)] . (3.19) 



By construction, the Hartree resummation does not change the first order result (</> 4 ) = 3 (4> 2 ) 2 ■ One may also check 
that it provides the correct next-to-leading term of order vA to (</) 4 )/A 2 , since the first, three-loop divergent diagram 
not of daisy type is only linearly divergent. 

Beyond the Hartree approximation the self-energy £(fc 2 ) receives k— dependent UV finite corrections which can be 
embodied into a finite k— dependent wave- function renormalization Z(k 2 ) by writing 

where ((f) 2 ) is now the exact expectation value and Z(k 2 ) is of course a function also of T and of the UV cutoff A, 
that is Z(k 2 ) = Z(k 2 ;T,K) 1 . We have from eqs. and ipHDJ) . 

1 _ fc 2 + l + £(fc 2 ) _ 1 t ±(k 2 ) 



Z(k 2 ) k 2 + l + 3((fj 2 ) k 2 +M 2 

where M 2 = 1 + 3 (4> 2 ) and £(fc 2 ) = £(fc 2 ) — 3 (4> 2 )- Notice that S(fc 2 ) can also be written as sum of Feynmann 
diagrams with dressed propagators (k 2 + M 2 )^ 1 and no tadpole insertion of any order. Thus, the lowest order 
contribution to S(fc 2 ) comes from the sunset diagrams with dressed propagators. Now, since the renormalized 4 
QFT is asymptotically free in the ultraviolet (that is the ultraviolet is controlled by the Gaussian fixed point), the 
exact 4> two-point function must behave as 1/fc 2 in the ultraviolet (even in the presence of the cutoff A, provided k 
remains much smaller than A). Thus we expect Z(k 2 ) to tend to unity for large k (and from above, since the two-loop 
sunset diagram enters in S with a minus sign), while (0 2 ) still diverges linearly with A. Since the only divergence in 
S(fc 2 ) comes from the sunset diagrams and is only logarithmic, Z(k 2 ;T,A) tends to 1 for all k in the limit A — > oo to 
any finite order of perturbation theory. A different behaviour in A would signal a nonperturbative effect. 

To know more about Z(k 2 ), in a nonperturbative way, we resort to Monte Carlo simulations in section IIII Dl 
In preparation, we need to re-express some relevant equilibrium quantities with the lattice regularization which 
corresponds to the discrete dynamics of section Hi CI 



C. Equilibrium on the lattice 



Thermal expectation values on the lattice are written as their continuum counterparts, eq. (|3.1I) . in terms of standard 
multiple integrals over phase-space degrees of freedom and the statistical weight exp(— H/T) where H is now the lattice 
Hamiltonian. As phase-space degrees of freedom we take 7r + (;c,0) = tt+(x) and (f>(x + aer + ,0) = <f)(x + a<r + ) with 
x G 2aZ 3 , or better x S 2o.Cn after the restriction to a finite volume (see section III Bp . As lattice Hamiltonian we 
take the total conserved energy of eqs. (|2.9|l and (|2.10() with D — 3, namely 



F[7r+,fl= ]T £+(x,0) = i^(2a) 3 [^ + (D0) 2 +0 2 h + 



¥1] 



i 



^(2a) 3 ir 2 + + (D0) 2 + \<5>^ + ±($« + a^ + ) 2 (l + ^) 



(3.21) 



where we used the first of eqs. Ij2.8|l to express <j)(x,a) as ^^(a:, 0) + cnr + (x,0). We see that the momentum field 
7r_|_ enters the lattice Hamiltonian in a way definitely more involved than in the continuum, eq. H2.4J) . or in the naive 
lattice regularization where only space is discretized and tt enters only the kinetic energy as in eq. I|2.6|) . However, 
H[n + , 4>] still depends on n + only ultra-locally, so that 7r + plays the role of a Gaussian auxiliary field which could be 
integrated away to produce an effective Hamiltonian for <f> alone with a non-polynomial local self-interaction. This 



Here A need not be a sharp cutoff as in the lowest perturbative order or in the Hartree approximation; our only request is that it is 
the UV cutoff in some spherically symmetric regularization procedure. For instance, one could use free-field propagators smeared in the 
ultraviolet by some smooth function of k/A dying at infinity faster than any power 
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self-interaction reproduces the standard continuum <fi A potential V[<fr] in the limit a — > 0. However, it is simpler not 
to integrate over tt + and keep working with the quartic Hamiltonian eci. (|3.21[) . 

We discuss first the Gaussian approximation to the thermal equilibrium on the lattice. This is valid for the lineariza- 
tion around classical solutions (zero field and the cnoidal) as well as for the selfconsistent Hartree approximation. 

Consider the lattice equation of motion (|2.12() and the first of the correspondence rules eq. I|2.8[l written in Fourier 
space, that is 

1 ° 
7f± lfc (t) = ±~[(t>k(t±a)-4> k (t)] , C(fea) = JJ cos %a (3.22) 

a 3 = 1 

Next we take the free-field limit of eq. I|2.12|l by linearizing it. Details are reported in Appendix [Bj where linearization 
is performed more generally over the uniform time-dependent solution of the discrete dynamics (the discrete cnoidal). 
Free-field dynamics is decoupled in Fourier space and harmonic also on the lattice, so that the standard virial theorem 
for harmonic oscillations (which holds true for phase-space averages as well as for time averages) yields, as in eq. (|H6f) 
specialized to D = 3; 

<'** |2 > = » 2 1 1 - tt|5 c M " i^' 2 > - [1+ i a ^i-W(M] s c ° , *" ,) (3 - 23) 

having used equipartition, that is (|7r+.fc| 2 ) = T/(l + \a 2 ), as follows from the quadratic part of the Hamiltonian 
ea. (|3.21|) . Go(k,a) in eq. (|3.23|) is the tree-level (Fourier transform of) the (j> two-point function on the lattice and 
can be seen to coincide with the result obtained by integrating 7T+ away within a quadratic approximation. Notice 
that (|0o | 2 ) = 1 as its continuum counterpart. The small a-dependent tree-level correction in (|7f +i fc| 2 ) = T/(l+ \a 2 ) 
w.r.t. the continuum (|7T;,| 2 ) = T is evidently a first effect of our specific lattice regularization. Similarly the finite 
lattice version of the tree-level tadpole reads 

,,2x 1 V- „ „ ^-ooT M , ml .. /-+^ 2 d 3 q 1 



Go(k,a) L =°°^[l + 0(a)]J 



L 3 ^ uv ' ' a 1 K "J _ /2 (2tt) 3 1-C 2 (q) 

ke{2i</L)c N J - 7r / 2 v ' KHJ (3.24) 

= 0.1741 ... - [1 + 0(a)] = 0.1108 . . .T A [1 + 0(A -1 )] 
a 

to be compared with the infinite- volume continuum expression eq. I|3.15fl . If these two expressions are identified, one 
obtain the tree-level relation A — 2.187 . . . A, which differs from that implied by (tt 2 ) [see section UlI A| . 

Another important issue about lattice effects concerns the equilibrium expectation value of </> 4 . We have already 
commented at the end of section lTl Cl on two different lattice regularizations of 4> 4 , the slightly nonlocal <j)± of eqs. I|2.8|l 
and the ultralocal </> 4 (cc) = [0(a;)] 4 . They do not yield the same values due to ultraviolet effects. At equilibrium this 
appears quite evident already at tree level: application of Wick theorem to the expectation value of the ultralocal 4> 4 
yields 

(0 4 )=3(0 2 ) 2 , (0 2 ) = i^Go(fe ) a) 

k 

as in the continuum, only with a sum over wave-vectors in the first Brillouin zone rather than a rotation invariant 
integral as in section Till Bl The same Wick rules on Gaussian integration plus translation invariance give instead for 

(0 4 ) = {($« +ai: + ) 2 ^) - [ } ' [[ 1 ' 



%a{\ + \a 2 ) (1 + ia 2 ) 2 



\2\ 



1 + ia 2 



T h 



8a 



2" 1 \*- 1 2 l 

2/f , I 



(1 + ia 2 ) 2 ' V 



2" 1 1 v 1 2" ) k 

A straightforward numerical integration yields the tree-level ratio 

/ {(j) 2 ) 2 = 1.154748... (3.25) 



for L = 12.8 and a = 0.0125. 



14 



When the interaction is turned on both (|-7f +j fe| 2 ) and (\4>k\ 2 ) gets modified. The relations between (|7r +j fe| 2 ) and the 
temperature T gets new corrections in terms of expectation values of dependent observables. A look at eq. I|3.21fl . 
with the knowledge that the Euclidean <fi 4 QFT is superrenormalizable (so that UV divergences can be fully assessed 
from perturbation theory as in the previous section), reveals that these corrections all vanish at least linearly as 
a — > 0, so that the continuum result (|7Tfc| 2 ) = T is recovered. At any rate, here we are only seeking a convenient 
parametrization for (|0fc| 2 ) in terms of (Itt+^I 2 ) and (</> 2 ). 

Again, we may resort to the linearization of the discrete equation of motion, this time over a uniform nonzero 
background <f>(t) , as detailed in Appendix [BJ First of all we learn how to obtain the Hartree resummation on the 
lattice: we replace the background <fi 2 (t) with the equilibrium expectation value (4> 2 ) in the dispersion relation eq. (|B4fl . 
obtaining 

1 + i a 2 (1 - u) 

cos[w fc a] = B{{(t> 2 ))C(ka) , B(u) = 2 (3.26) 

[l + ia 2 (1 + u)] 

and invoke the virial theorem (the Hartree approximation is Gaussian, that is describes harmonic oscillations) to write 

nl |2v _ fl2 (r^+.fcl 2 ) a 2 (l^+,fc| 2 ) n97] 

; 1 - (1 - \ a 4 )cosW ~ 1 - [S((0 2 )) C(ka)} 2 [ ' 

We have neglected the term a 4 cos 2 ojk in the last step since it is of order a 4 for all k. Beyond the Hartree approxima- 
tion, we promote eq. (|3.27[1 to an exact relation by introducing the k— dependent lattice wave- function renormalization 
Zk = Z k (T,a) 

n 7 |2\ a 2 Z k (|7T+,fc| 2 ) . . 

(l0fc| >= l-[B(m)C(ka)} 2 (3 - 28) 

Zk(T, a) is by construction as much as possible free from lattice effects (such as the deformation of free-field propagator 
at large |fc| due to the replacement of continuum derivatives with finite differences) and, once averaged over discrete 
directions to Zk(T,a), it should provide a rather accurate representation of its continuum counterpart Z(k 2 ;T,A), 
if a universal relation between A and A = w/(2a) could be determined. As we have seen above, in the case of UV 
divergent quantities like (tt 2 ) and (</> 2 ) this is not possible, but should be possible for renormalized or UV finite 
observables. Perturbation theory suggests this to be the case for Z(k 2 ;T,A). In the next section we shall follow a 
different, nonperturbative strategy, based on the simple observation that the agreement of lattice calculations with 
their continuum A— cutoffed counterparts is guaranteed for all observables in a different limit, that is A — > oo first, 
at fixed large A, rather than A, A — > oo with A = 0(A). For k— dependent quantities this means agreement for 
1*1 < A « A. 



D. Monte Carlo simulations 



Monte Carlo simulations allow to compute thermodynamical averages beyond perturbation theory and Hartree 
resummation. We performed Metropolis simulations for the </> 4 model discretized on the lattice as in section Hi CI 

Consider the lattice Hamiltonian H[tt + , 4>] of eq. (|3.21|) . In this context it is more convenient to revert to the notation 
with the single field <fi, using both 4>(x,0) and <p(x,a). We may drop the distinction between the two time slices 
by considering the degrees of freedom attached to a cell-centered cubic lattice. There are therefore the points of 
2aZ 3 , labeled by x = 2an, n having integer coordinates, and the points of 2a(Z + 1/2) 3 labeled x + acr + , where 
<r+ = (1,1,1). The Hamiltonian can then be written as 

# = 2a(2 + a 2 ) ^ [cj) 2 {x) + <j) 2 (x + acr + )] - a ^ ^ (j){x) (j){x + aer) 1 - \a 2 (j>(x) (j>(x + acr) . (3.29) 

This Hamiltonian contains local terms and couplings between nearest-neighbors along the main diagonals. The quartic 
term contribute to these couplings and there is no quartic local term. However, to perform Monte Carlo simulations is 
more convenient to add and subtract a quartic local piece. In this way one can change the local integration variables 
— oo < 4>{x) < +oo in Eq. (|3.1|l to a new variable C(x) ranging from zero to one as in [29l |. 
We therefore split our Hamiltonian eq. (|3.29|l as 



H = H 1 +H 2 , 
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# l = E {hi[4>{x)] + h l [(i){x + acT + )]} , h x [4>] = 4 a (2 + a 2 ) cj) 2 + 2 a 3 (/> 4 , 



xG2aZ 3 

,3 



H * = - a E E^^ + ^ + t E E^)^ 2 ^^)-^)] ■ ( 3 - 3 °) 



4 

x(E2aZ 3 <t x<E2aZ 3 <r 



The functional integral in Eq. l|3.1|l can be written as, 



(9) 



Jo • • • Jo e[c(.)] e-^ra 



where the new field C is defined as 



(3.31) 



dxe-^M 

J — OO 

f +QO dxe-P 

J — OO 



C(^) == 7+^ JTTT ■ (3-32) 



In this way the classical </> 4 model is mapped into a classical continuous Ising-like model where the dynamical variables 
C(x) run between zero and one. The Ising-like Hamiltonian is given by if 2 [Eq. Ij3.30|l ]. 

Following the Metropolis method we generate a sequence of configurations for the whole cell-centered cubic lattice, 
as follows. We start by choosing random values at each point of the lattice for the variables C(x). Then, we choose 
one point in the lattice x at random and consider a new value C for it picked at random between zero and one. We 
then have to go back from the variables C(x) to the variables 4>(x) in order to compute the energy of the old and 
new configurations. Notice that the inverse function <j> = 4>(C) is unique since dC/d(f> > as one sees from Eq. (|3.32|l . 
Their energy change is given by 

A = # 2 (new) - #2 (old) =2a [<j)' - cf>(x)} <p 2 {x + acr) - 4 [(j)' 2 + <p 2 {x)} j 

where we used Ea. (|3~3UJl and <p' = 4>{C") 

We now follow the standard Metropolis procedure. That is, we compare e _/3 A with a random number r with 
< r < 1. If r < e~P A we pick the new configuration. Otherwise, for r > A we keep the old one. We repeat this 
process many times producing in this way a sequence of configurations on which we can compute the values of any 
observable 6. 

It can be shown (2^| that the average of these values of <d converges to the thermodynamical phase average in the 
limit when the number of iterations AT goes to infinity The expectation value is then given by 

M 

< e > = ^E «- 

n— 1 

In the continuous limit a — > we see from eq. (|3.30|l that the quartic terms become negligible and hence all quantities 
will be functions of f3 a — a/T. In addition, gradients of the fields scale as l/a. 

In fig. ^ we plot (0 2 )j ((V0) 2 ) and the ratio (</> 4 )/ K^ 2 )] 2 as functions oiT/a from Monte Carlo simulations. We see 
that ((j) 2 } and ((V0) 2 ) follow the linear behaviour predicted by low temperature perturbation theory eq. (|3.15() while 
the ratio (4> 4 ) / [(<fi 2 )] 2 takes the value 3 predicted by Wick's theorem in low temperature perturbation theory eq. (|3.19|) 
or in Hartree approximation. In the Monte Carlo simulations we used in the lattice the ultralocal definition for (jr. 
For the definition w e find from the Monte Carlo simulations (4>\.) /[(4> 2 )} 2 — 1.1 ... for a = 0.05 and L = 3.2 in 
agreement with ea. (|3.25|) . 

As a matter of fact, by measuring also the two-point function {<f>(x)<f>(x')) we verified that the Hartree approximation 
is indeed quite accurate for all wavelengths, except for the largest ones, and for a wide range of temperatures. We 
recall that this fact depends on the classical statistical mechanics interpretation of the ultraviolet divergences, which 
are physical and cannot be eliminated by counter-terms as in quantum field theory. Thermal equilibrium at nonzero 
temperature really makes sense for the classical </> 4 theory in three space dimensions only with a UV cutoff present, that 
is on a lattice. As a — > 0, the only way to keep things finite is to decrease T — > 0, making the Hartree approximation, 
which resumes all T/a corrections, more and more accurate. 
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FIG. 1: {4> 2 ), ((V(j>) 2 ) and the ratio {<f> 4 } / (4> 2 ) 2 as a function of T/a in thermal equilibrium from Monte Carlo simulations. 



IV. DYNAMICS OF THERM ALIZATION 

We present here the exact evolution of the lattice 4 theory defined in sec. Ill CI and the development of the 
ultraviolet cascade. We considered three system sizes, L = 6.4, 12.8 and 25.6 and several lattice spacings, a = 
0.1, 0.064, 0.05, 0.025, 0.0125 and 0.00625. The largest system contains 1024 3 lattice sites. We worked with the 
largest volume 25. 6 2 only for a = 0.05 and E/V = 89.5, verifying that the change from L = 12.8 to L = 25.6 does not 
change any observable in a significant way. Therefore L — 12.8 has been our preferred choice. 



Initial Conditions 



We used a large variety of initial conditions in our calculations. In these studies the initial power is concentrated 



in the infrared; that is, |0fc| 2 (O) and | Trjt | 2 ( ) are significant only for wavenumbers well below the cutoff A 
particular, we considered superpositions of infrared plane waves, when the initial fields have the form 



A" 



(x, 0) = 0o + A 2_, c i cos(fci • x + 2 ir-fi) , w(x, 0) = B \_. d% cos(fc.; ■ x + 2 ir Sj) 



i=l 



as well as superpositions of localized wave packets of the form 

<fi(x, 0) = 0o + A a w(x — Xi) , it(x, 0) = B di w(x — x[) 



In 



(4.1) 



(4.2) 



where 



w(x) 



w (k [x - L n] ) 



and wq(x) is either the Gaussian, Wo(x) = e x , or the Lorentzian, wq{x) = (1 + x 2 ) 1 . The sum over n in the last 
equation is needed by p.b.c, but in practice, with our choice L ~ 10, only few terms in the sum are needed. These 
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fields have support throughout Fourier space, but peaked as Gaussians or simple exponentials at low wavenumbers 
h < h 

' L rsj "-max- 

In eqs. (|4.1|l and 14.2|l . (f>o is a uniform background, or homogeneous condensate, while the wave- vectors ki = 2irrii/L 
in eq. 14.1fl have nonzero modulus |fc| < k max < 30ir/L — 15n/(Na) <C n/(2a). The number K in ea. (|4.1|) is chosen 
within the range 10 — 100, proportionally to fc max , and the specific ki in eq. 14.11) or the Xi and Xi in in eq. I|4.2[) are 
chosen at random. The phases 7, , Si in eq. (|4.1(l and the relative amplitudes Cj , di in both cases are chosen at random 
in the interval [0, 1]. Finally, for any given realization of these random numbers, the background 4>q and the overall 
amplitudes A and B are constrained in such a way that the energy density E/V takes a given, predefined value. 
Typically, we further restricted the remaining freedom by choosing either B — or B — A, so that the 'macroscopic' 
part of the initial state is entirely characterized by the values of E/V, of the condensate 4>q and of fc max - In principle, 
one should then regard the randomly chosen numbers as 'microscopic', or fine-grained properties providing the initial 
entropy; one should then average all measured observables over these random choices, as described in section Hi DI to 
reduce the fluctuations in the measures. Actually, we performed such average (over 20 to 40 initial conditions) only 
rarely, since fluctuations were kept under control already by the 3D space average, by the sliding time average and/or 
by the average over discrete directions of Appendix lAl 

It should be noticed that all our initial configurations have vanishing total momentum J d 3 x ir(x) V0(a;), which is 
a conserved quantity for p.b.c. 



B. Time evolution of local observables 



Here we show the time evolution of the basic local observables of eq. H2.16[) averaged in space and time according 
to section Til DI As mentioned above, we considered many different initial conditions, with several values of the energy 
density ranging from E/V = 0.1 up to E/V = 5000. 

Figs. |5]to[Sl display the local observables 7r 2 (£), (V</>) 2 (i), <fi 2 (t) and 4>%{t), as functions of time for different values 
of the relevant parameters. When an initial zero-mode condensate is present [ (f> in Eas. l|4.1|l and l|4. 2f> ], we quantify 
its weight with the ratio Eq / E between the energy built solely from the zero- mode and the total energy. 



log 7T 2 




4.5 - 



L — 6.4 
a = 0.1 
E/V = 1000 




4 5 6 7 8 9 10 11 12 

logt 



FIG. 2: Log-log plot of the time evolution of local observables with sliding time average. Initial conditions are localized wave 
packets of Lorentzian shape without any zero-mode condensate. No average is performed over initial packet amplitudes or 
positions. 
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FIG. 3: Log-log plot of the time evolution of local observables. The time averaging interval is quite small, r = 2, and the initial 
conditions are plane waves with a large zero-mode condensate. No average is performed over initial amplitudes or phases. 
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FIG. 4: Log-log plot of the time evolution of local observables. The time averaging interval is r = 40 and the initial conditions 
are plane waves with a large zero-mode condensate. No average is performed over initial amplitudes or phases. 
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FIG. 5: Log-log plot of the time evolution of local observables. Initial conditions are plane waves without zero-mode condensate. 
The time averaging interval is r = 200 and an extra average is performed over 40 random choices of initial amplitudes and 
phases. This explains why the curves in this picture are smoother than those in fig. [3] and 0] 

These plots clearly reveal the main qualitative feature of the thermalization process in 3 + 1 dimensions. Initially 
(V0) 2 (£) is small reflecting the fact that the initial conditions determine a power spectrum localized at wave- vectors 
with k -C 1/a. The mode mixing entailed by the interaction is transferring power to larger wave- vectors, thus 
effectively transferring energy from the interaction and 7r 2 terms, which decrease, to the spatial gradient term which 
increases. As is clear from these figures, Tr 2 (t) and (V</>) 2 (i) tend to a limit for late time. (j> 2 (t) and 4>%{t) show a clear 
limiting behaviour only for E/V large enough. 

The steady growth of (V0) 2 (t) at the expense of the interaction term 4>\{t), the kinetic energy 7r 2 (i) and the mass 
term (f> 2 (t) shows that the basic mechanism leading towards thermalization is a rather uniform flow of energy towards 
larger k modes, namely the ultraviolet cascade. We can associate a time scale to to the onset of this ultraviolet cascade 
by looking at the point when (V</>) 2 (i) overtakes 4>%{t), if such point exists. For E/V too small, implying small field 
amplitudes and very small initial 0+(O), (V0) 2 (t) might be larger than <j)\if) from the beginning. In this case we can 
take the point where the rise of (V0) 2 (i) is the steepest. ^From figs. |21to|51 as well as from the rest of our data, we 
see that to does depend on E/V and on the condensate, but stays within the range 100 < t < 400 for E/V in the 
range from 10 to 5000, with to decreasing as E/V increases at fixed value of the condensate. 

We are interested in somehow large energy densities E/V which is the relevant regime both for the early universe 
and the ultrarelativistic heavy ion collisions. However, it is physically interesting to also study the low energy density 
regime. We find that thermalization happens for small values of E/V . No threshold to a non-ergodic behaviour was 
found. However, the thermalization dynamics slows down dramatically when E/V is reduced well below a value ~ 10 

and t grows substantially. We depict in fig. (V</>) 2 (T), <j) 2 (t), lj?{t) and ^(i) for E/V = 0.1, a = 0.1, L = 6.4. 
We see that the ultraviolet cascade only starts in this case for much later times \nt ~ 11, t ~ 50000. 

Fig. [^displays log | </>(£) | as a function of the logarithm of time for E/V — 1000, 10 and E/V = 0.1 for L = 6.4 and 
a = 0.1, respectively. We see that the relaxation of (j) towards its thermal equilibrium value (</> = 0) is different from 
the other physical quantities previously discussed. This is due to the fact that the vanishing of <fi is connected to a 
symmetry of the model. We find that <fr(t) relaxes as ~ 1/t for In t > Into — 6, t > to — 500 in the universal stage. 
This can be seen from fig. |SJ 
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FIG. 6: log \<j>(t)\ as a function of the logarithm of the time t for E/V = 1000 and E/V = 10 with L = 6.4 and a = 0.1 . No 
average is performed over initial packet amplitudes or positions. Time averaging is performed here as by eas. 12.171 - 12.181 . 

We plot in fig. 08 (</*+) vs. (4> 2 ) in thermal equilibrium (obtained from Monte Carlo simulations) as well as <f)\{t) vs. 
4> 2 (t) from time averages. We see that the curve obtained from the lattice Monte Carlo calculations is in agreement 
with the ones from time averages. Since both (j>%(t) and (j) 2 {t) vary with time, the agreement with the thermal curve 
implies that, at least for these observables, we are in a situation of effective thermalization with a time depending 
temperature. As in ref. the effective temperature decreases with time. The small disagreement (less than 10%) 
between the thermal equilibrium curve (from Monte Carlo) and the time averages in fig. |S| comes from the low k 
modes contribution. As it will be discussed below, infrared modes are much slower to thermalize than the modes with 
k 2 > cf> 2 (t). Actually, the disagreement decreases with increasing time [decreasing <fA_(t) and <j> 2 (t)] showing that the 
1R modes are getting thermalized equilibrating with the modes with k 2 > 4> 2 (t). 



C. Virialization and Equation of State 



We depict in fig. [5] the quantity, 



A(t) = 



(^ 2 )(t)-((V^) 2 )(t)-<0 2 )(t)-<0l.)(t) 
«* 2 >W 



(4.3) 



This quantity vanishes when the virial theorem is fulfilled [see ea. <f3.6f) ]. It turns out to be negative for finite times 
and nonzero a. We see from fig. I^Jthat |A(i)| starts to decrease at times earlier than to. Therefore, the model starts 
to virialize before it starts to thermalize. |A(t)| keeps decreasing with time and tends to a nonzero value which is 
of the order 0(a 2 ) for t — > oo. This is to be expected since ea. H3.6fl only holds in the continuum limit and receives 
corrections in the lattice 0(a 2 ) as shown in ea. H2.13fl . 

We computed the pressure as a function of time from eq. (f3.10l) as 



P= 2 



(V0) 2 



Notice that we are not using the virial theorem. 
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FIG. 7: Log-log plot of the time evolution of local observables with sliding time average and very small energy density. Initial 
conditions are of the same type as in fig.|U No average is performed over initial packet amplitudes or positions. Time averaging 
is performed here as by eqs. 12. 171 - 12. 181 . 




FIG. 8: (4>+) vs. (4> 2 ) in thermal equilibrium (obtained from Monte Carlo simulations) and 4>\(t) vs. (f> 2 (t) from time averages. 
Both 4>\(t) and (j) 2 (t) decrease for increasing time. 
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FIG. 9: The normalized l.h.s. of the virial theorem A(t) vs. the logarithm of the time t for E/V = p = 1, 10, 100 and 1000 
for a = 0.1 and L = 6.4 [see eq. 14.3ll l. No average is performed over initial packet amplitudes or positions. Time averaging is 
performed here as by eas. l|2.17|l - l2.18^ . 
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We depict in fig. EHf/p as a function of time. We see that 



for the whole range of p considered. This inequality is in agreement with ea. l|3.11|) . That is, fig. 1101 shows that the 
equation of state approaches approximately the radiation-dominated equation of state unless p is too small. 

Notice from figs. 12151 and fig. that the approach to the radiation-dominated equation of state parallels the 
decrease of <fA since both are governed by the UV cascade. 
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FIG. 10: The equation of state p/p vs. the logarithm of the time t for p = 100 and 10 with a — 0.05 and for p = 5 and 0.1 
with and a = 0.1. We have in all four cases L = 6.4. No average is performed over initial packet amplitudes or positions. Time 
averaging is performed here as by eas. 12.17H - 12.18ll . 

We see that the virial as well as the equation of state reach stationary values despite thermalization is not achieved. 
As already noticed in ref. the equation of state follows by taking time average over the period for the cnoidal 
solution. Moreover, as we show in the Appendix, the virial theorem is exactly fulfilled by the cnoidal solution averaging 
over one period. These type of phenomena has been recently highlighted in ref. |24j | . 



D. Evolution of power spectra and UV cascade 



We now turn our attention to the study of correlation functions. 

According to the Fourier transform relationship cq. 1)2. 21|) between the power spectrum |</>fc| 2 (i) and the equal-time 
correlation function <jxj){x,t), there are two approaches to the numerical evolution of such quantities (we specialize 
now on 7r but the discussions applies equally well to </>). We extract the field 7r(x, t) from the lattice fields F(n, s) and 
G(n, s) [see first line in eq. (|2.8J) and eq. H2.14J) ]. Fourier-transform it to TTk(t) and then perform all needed averages 
on |7Tfc(<)| 2 . Or we directly compute averages of the correlations of F(n,s) and G(n,s) and extract from them the 
correlations WW(x, t). We found that both methods yield the same results. 

Moreover, when using the approach with growing time averages as in ea. l|2.18|l with a unique initial condition, one 
realizes that the simply time-averaged correlation 

1 r* 

- / dt' 7r(sc, t') ir(x', t') (4.4) 

T Jt-T 
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very soon (in the logarithm of time) becomes translation invariant, that is a function only on the distance \x — x'\, 
making the time consuming space average unnecessary. 

Figs. II II and IT2l show the average power |7r/c| 2 (£), multiplied by the spherical measure Airk 2 , for five values of the 
lattice spacing and one given choice of all other parameters. The region where the power is significant is spreading 
towards the UV cutoff. 
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A = 62.832 



FIG. 11: The power Ank \TTk\ 2 (t) vs. k = |fe| at 20 different times ranging, in an approximately exponential way, from t = 
to t — 12.75 (left) and from t = 12.75 to t = 2948 (right). The time averaging interval is r = 2 and the initial conditions are 
infrared random plane waves, as apparent from the IR peaks at early times. The front of the UV cascade arrives close to the 
cutoff A for the latest times. 



The chosen initial conditions eas. (|4.1|) - (|4.2|) are such that the power is concentrated in long wavelength modes 
with k well below the ultraviolet cutoff A = 7r/(2a). Therefore, | tta; | 2 (0) is concentrated on small k. During the time 
evolution the non-linearity gradually transfers energy off to higher k— modes leading to the ultraviolet cascade as 
discussed above. 

It is important to observe the effects of the finite UV cutoff A. At the given value E/V — 89.5 of the energy density, 
when A = 62.832 (Fig. 111(1 . which corresponds to a lattice with 256 3 sites, the shape as a function of k and the time 
evolution of the power spectrum are still markedly distorted by the presence of the cutoff, although much less than 
for the smaller values A = 15.708 and A = 31.416 (Fig. [TJJ). Only when doubling the cutoff from 125.664 to 251.328 
(corresponding resp. to 512 3 and 1024 3 lattices) the effect of the cutoff in |7Tfc| 2 (i) does not appear significant. 

That is, the cascades shown in the lower panels of fig. 1 121 are cut-off independent since they are evolving well below 
the cut-off A. 

Plots analogous to figs. 1111 and 1121 can be drawn for the field power k 2 \<j>k\ 2 {t)- However, since unlike k 2 |frfc| 2 (£), 

k 2 |<^fe| 2 (i) never grows with k [at thermal equilibrium and in the continuum |^fc| 2 (t) goes like k~ 2 , see section Till Bj . 

the UV cascade for k 2 \(j>k\ 2 (t) is not as evident as for k 2 \TTk\ 2 (t). The time evolution of both powers can be better 
appreciated in figs. 1131 and !14l One can see a rather complicated behaviour in the strongly interacting infrared, while 
the rest of the modes exhibit a rather orderly evolution which indicate a weakly interacting dynamics. 

One sees from figs. ED an d El that m0( i es with low k decrease in amplitude with time except for k = and the 
first two non-zero modes. Modes with larger k grow in amplitude monotonically with time showing the existence of 
the smooth UV cascade. Since total energy is conserved the larger k modes grow in amplitude at the expense of the 
lower k modes. Notice that only a restricted number of low k modes have a significant initial amplitude. They feed 
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FIG. 12: The power 4tt k 2 \Tv k \ 2 (t) vs. k as in the right of Fig. E] but at values of the UV cutoff A scaled by 1/4, 1/2, 2 and 
4. When A = 251.328 the last time t — 2978 is missing. The cascades shown in the lower panels are cut-off independent since 
they are evolving well below the cut-off A. 




1 2 3 4 5 6 7 8 

log* 



FIG. 13: Log-log plot of |<j!>fc| 2 (t) vs. time. Initial conditions are as in fig. 1111 The four thicker lines correspond to the modes 
initially filled (after averaging over directions). The zero-mode (dashed line) was not filled at t = 0. Notice the peculiar 
behaviour of the lowest fc-modes compared with the rest of the fc-modes. 
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FIG. 14: As in fig. 1131 but for |7Tfc| 2 (t) vs. time. Notice the peculiar behaviour of the lowest fc-modes compared with the rest 
of the fc-modes. 



the growth of a large number of modes with larger k which start with very small initial amplitudes. Later, for t > to 
all IR modes decrease with time. 

As said before, the k = mode and the first two non-zero modes in figs. 1131 and IT~i1 behave differently to the rest of 
the modes. They start by growing with time and they stay larger than all the other modes for a while. Actually, 

as we shall see below in sec. IIV II . the modes with < k < \j~<j?{t) behave differently to the rest keeping a significant 

coupling among them while modes with k > \J~<fi 2 {t) exhibit weak nonlinearities for late times in the lattice model. 
More precisely, these latter modes obeys the Hartree approximation and exhibit effective equilibration much earlier 
than the infrared modes. 



A single parameter that efficiently measures the UV cascade for |7Tfc| 2 (i) is the average wavenumber k(t), defined 
[in continuum notation and recalling the sum-rules eq. I|2.20p ] as 



7T»(i) 



d 3 k 



\k\W(t). 



(4.5) 



However, it will be more appropriate to exclude from the integration in ea. l|4.5|l the infrared modes. More precisely, 
since for the evolution time considered, a large portion of energy lingers on the IR modes filled at t — (see fig. IT^ . 
one should restrict the averaging over the modes not filled at t = and define more properly, after averaging over 
discrete directions, 



*(*) = p 



k 2 dk 

2 IT 2 



k\n k \ 2 (t), P 



k 2 dk-. 



2tt 2 



(4.6) 



where ko is the smallest unfilled radial wavenumber larger than all wavenumbers of the modes filled at t — 0. This 
procedure applies directly when the initial conditions are superposition of IR plane waves. In case of superposition of 
localized wave packets, k may be defined as the value of k at which |7ffc| 2 (0) drops below some fixed small value. 

We plot k(t) in fig. E3for the same context of figs. 1111 and IT21 The cutoff effects at larger values of a and the 
saturation in the continuum limit a — > are quite evident. Still, in spite of the fact that the UV cascade has time to 
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fully develop, as evident from figs. lllland lT2l fig. 1151 shows that the growth of k(t) with time is not a simple power, 
not even for the latest times considered. However, 

k(t) ~ k t 1 / 3 (4.7) 

provides a rough estimate. 



3.5 



o a 


= 0.1, logfccq = 2.714 


a 


= 0.05, logfc cq = 3.407 


v a 


= 0.025, logfc cq = 4.100 


— * — a 


= 0.0125, logfc cq = 4.793 


a 


= 0.00625, logfccq = 5.486 



g>2.5 



1.5 - 




FIG. 15: Log- log plot of the time evolution of average wavenumber k(t), for five values of the lattice spacing 2a, with initial 
conditions as in figs. llll and [T2l 



Let us now turn our attention to the behavior in k of |7Tfc| 2 (i) during the UV cascade. From all our data, it is clear 
that |7Tfc| 2 (t) dies exponentially fast for k > k(t) as long as k(t) <C A [sec fig. [1^1 as an example]. This exponential 
behavior is clearly time-dependent. For k not too small nor too close to k(t), |7Tfc| 2 (t) exhibits a decreasing power- like 
behaviour, as can be seen for instance in fig. 1171 where the log-log plot of 47r k 2 |7Tfc| 2 (i)) vs. k is shown. That is, 
|frfc| 2 (t) ~ k~ a at different times, where the values of —a are obtained subtracting 2 from the numbers indicated in 
fig- El ol decreases monotonically from 1.12 at t = 812 to 0.17 at t = 48385, when the forefront of the cascade is 
reaching the cutoff A. 

In figs.^]and^]we give some graphic examples of the UV cascade in the presence of an initial zero- mode condensate 
[4> in Eas. l|4.1|> and l|4.2|) ]. We notice the evidence for the parametric resonance, whose location in A;— space well agrees 
with the analytic prediction from the solution of the Lame equation for the linearized model [see ea. l|B7|) ]. It is quite 
evident however, that such parametric resonance plays practically no role in the UV cascade, as already noticed in 
ref.^^|. We also stress once more the importance of UV cutoff effects when a is doubled from a = 0.0125 to a = 0.025 
(corresponding to a reduction of the lattice from 512 3 to 256 3 sites). This is particularly relevant since the zero-mode 
(and few lower k modes) remain quite large for all times considered (see fig. I19[): according to the scenario of ref . [Isj . 
the regime should then remain that of driven turbulence until the lattice effect become dominant; thus, for the energy 
density of figs. 1181 no cutoff-independent regime of free turbulence is observable in our evolution when a = 0.025 
and L = 12.8. (Notice that our lattices reach a size 4 3 times larger than those in ref. 18]). 

It would be indeed very interesting to find this universal UV cascade in quantum th eory for large occupation numbers 
and small couplings. Most of the works consider other regimes 

llilialia- In ref. [16| the small coupling regime m 
quantum theory is investigated including the leading 1/N corrections. However, thermalization is not reached in ref. 
[ll| since the times considered are not long enough. The classical evolution is compared with the quantum evolution 
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FIG. 16: log [|-7Tfc | 2 (t)] vs. k for several values of time. All relevant parameters are at the indicated values. Initial conditions 
were random plane waves. |7t"fc| 2 (t) starts decreasing with k as ~ k~ a and later it dies exponentially for A > k > k(t). 
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FIG. 17: log[4 7r fe 2 1 TTfe | 2 vs. logfe in a fixed wavenumber window for several values of time. All relevant parameters are at 
the indicated values. 



29 



to first order in 1/N for the <fi 4 model in 3 + 1 dimensions in ref. [23. It is stated there that the classical evolution 
is a good approximation to the quantum evolution for non-asymptotic times therefore supporting the relevance for 
quantum field theory of the classical dynamics studied here. 
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FIG. 18: The power 4ir k \nk\ 2 (t) vs. k — \k\ at several times, E/V = 89.5 and two lattice spacings. The initial conditions 
are infrared random plane waves with a zero-mode condensate dominating the energy. Notice the peaks due to parametric 
resonance. 



E. Thermalization of the power spectra 



We present here the results of a long simulation yielding to a very good approximation thermalized powers |7r/t| 2 (i) 

and \(j>k\ 2 {t), on a lattice cube C N with N = 100, a = 0.064 and L = 12.8. The energy density is E/V = 569.5, still 
relatively small compared to A 3 ~ 14785. The initial conditions are random plane waves as in eq. I|4.1|l . with the 
initial powers of cf> and ir exactly null for k > 3.436. The rather large time averaging interval t — 200 is still negligible 
compared to the length t = 915494 of the evolution. To further reduce fluctuations, in this case we also performed an 
average over 32 different random initial mode amplitudes and phases. 

In Fig. [3U|we plot the basic one-point observables and the average wavenumber k(t) as function of time. The 
system clearly reaches a time-independent stage for \nt > 12.5. In Fig. |2]we plot the power spectra |7Tfc| 2 (i), times 
the spherical volume Airk 2 , vs. the wavenumber k at several times: the UV cascade up to the cutoff is evident. We 

plot in log-log scale both |7Tfc| 2 (i) and |0fc| 2 (<) vs. k in Fig. 1221 and vs. time in Fig. [23 Again, the system shows an 
evident limiting behaviour: the large initial peaks in the IR modes have almost completely disappeared for logt ~ 12.5 
and at the latest times logt ~ 14 the power spectrum |7Tfc| 2 (i) is flat, except for fluctuations and a small remnant 
of the infrared peaks exactly at k — (see top of fig. I24|) . From the height of |7Tfe| 2 (t) in the plateau we read the 
temperature T = 1.2, very close to E/N 3 — 1.19432 . . .. The difference is almost completely accounted for by the last 
value of ^ = 10.264 . . . [see eq. (I3~H|) ]. 

To ascertain the thermalization of |<^fc| 2 (t) is very convenient to consider the time-dependent analog Zk(t) of the 
equilibrium wave-function renormalization Z\~ [see ea. ()3.28[) ] replacing there the equilibrium quantities by their time- 
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FIG. 20: Log-log plot of the time evolution of space-averaged local observables and of the average wavenumber. The time 
averaging interval is r = 200 and the initial conditions are infrared random plane waves. This plot is analogous to fig. H but 
with the time evolution lasting for significatively later times. 
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FIG. 21: 47rfc 2 |7Tfc| 2 (t) vs. k at the times indicated. All parameters are as in Fig. 1201 





FIG. 22: Log-log plot of |S"fe| 2 (t) and |0fe| 2 (t) vs. k. Data points are explicitly indicated. All parameters are as in Fig. 1201 
Comparison with fig. 1201 shows that modes with k 2 < <j) 2 {i) thermalize much slower than those with k 2 > (j> 2 {t). 
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dependent out-of-equilibrium counterparts, that is 



\M 2 (t) = 



a 2 Z k Jt)\* k \2(t) 
l-[B(^(t))C(ka)f 



(4.8) 



This relation implies that the mode k has indeed virialized in the interacting theory. 

In the bottom of fig. [21 we plot the (average over discrete directions of) the equilibrium and Zk(t) for several 
late times. It is evident that Zk(t) approaches its equilibrium value Z\. earlier and much more closely than |7rj,| 2 (i). In 

other words, when both |7rfc| 2 (i) and |<?i>fc| 2 (i) are still out of equilibrium over a wide range of low wavenumbers, Zk(t) 
is practically already at equilibrium except for a very narrow range of very small wavenumbers. One can see that this 
range is approximately given by k 2 < 4> 2 (t). These infrared modes are the last to effectively thermalize. [Recall figs. 
I13ll4l and the discussion in sec. IIV Dl about the peculiar behaviour of these low A;- modes]. 

In conclusion, the field does thermalize completely in the lattice with a time scale of the order 10 6 . For such 
times |7ffc| 2 (t) is k independent as seen in fig. [521 Contrary to the 1 + 1 dimensional case[fj, the infrared modes with 
k 2 < (j> 2 (t) thermalize here the last. 





F. Equilibration in the cascade 

The effectiveness of the time-dependent wavefunction renormalization Zf.(t) in the study of the equilibration process 
for very late times, suggests to use it also at earlier times, during the universal UV cascade. As a matter of fact, 
we may regard eq. (|4.8() just as a definition of Zk(t), which is possible for any time but might very well turn out to 
yield a function of k quite different from its equilibrium counterpart Z k . What we find, however, is something very 
close to equilibrium already at times of order to ~ 500, which are those proper of the universal cascade. To have a 
cutoff-independent window wide enough, we consider here a lattice spacing a — 0.0125, much smaller that that of the 
the previous section, and a smaller energy density E/V = 89.5. Notice that this would correspond to an equilibrium 
temperature of order E/N — 0.001398 . . ., almost negligible as compared to that of the previous section. 
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FIG. 24: Comparison of log \nk\ 2 (t) vs. Zk(t). Diamonds in the bottom represent equilibrium values. All parameters are as in 
Fig. 1201 We see that Zk{t) equilibrates much earlier than |7ffc| 2 (t). 



In Fig. [23 on the left column, we plot Zk(t), the direction averaged Zk(t), for several times ranging from t = 135 to 
t = 2978 and for 0.695 ...< k< 9.559 ... in the top and 4 .644 ...<k< 28.218 ... in the bottom. On the left column, 

we plot Zk(t) vs. the scaled wavenumber u = kj \j 4> 2 {t). The good collapse in the bottom left shows that Zk(t) is 
a almost function only of u in the range 3.430 ... < u < 20.843 . . .. The collapse for 0.513 < u < 2.5 is definitely 
worse, in agreement with the slower equilibration of the infrared modes. Most remarkably, by direct comparison with 
fig. 1241 we observe that Zj~(t) is very similar, qualitatively and quantitatively, to its equilibrium counterpart at a 
temperature roughtly a thousand times larger than the temperature the system will eventually reach when t — > oo. 
One could say that equilibration in the bulk of cascade has taken place with an effective temperature of order one. 
The precise time dependence of this effective temperature as well as a check on its universality (that is to say unicity 
when the studied observables are changed) require a much more elaborate analysis which is beyond the scope of the 
present work. Here we only stress that, as a function of u, Zk(t) resembles closely its equilibrium counterpart and in 
particular it decreases to unity rather fast as u grows, that is when k 2 3> <j) 2 (t). 



G. Effective frequency and particle number 

The motion of each Fourier mode <pk of the field is characterized by fast oscillations into an envelope which varies 
relatively slowly in time. In the linear approximation, the scale of fast oscillation is fixed by the frequency of the 
free massive dispersion relation, that is ujk = k 2 + 1 in the continuum. Of course on our staggered lattice one should 
consider the lattice dispersion relation, which differs from the continuum form by power corrections in a 2 k 2 . In 
Appendix [B] we show that the exact free dispersion relation on our lattice reads 

1 3 

cos(wfca) = j— TT cosfcj-a (4.9) 

1 ~\~ ~^G> 

2 J = l 

which correctly reproduces k 2 + 1 in the limit a — > at fixed k. 
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FIG. 25: Scaling properties of the k— dependent wavefunction renorrnalization for two low-lying ranges of the radial wavenumber 
k. 



On the continuum the modes <pk obey the Fourier transform of the exact field equation (|2.3I) , namely 

d 2 



dt 2 



bk = -ujk 4>k - (4> )fc 



(4.10) 



where {<fi 3 )k is the Fourier transform of (j> 3 (x). 

Now to disentangle fast and slow motions one can use time averages over intervals large compared to the periods 
of the fast oscillations but small compared to the time scale of the slow motion. This is exactly what we did above, 
revealing the smooth UV cascade. To give a quantitative description of this scenario we can use the following simple 
argument. Consider the quantity 



-k^k 



If the motion of each mode were strictly periodic, the time average over several periods of Ik would vanish. But this 
is also true if the fast and slow motions are indeed separable, at least after some coarse graining like averaging over 
discrete directions. Then we may proceed as in the standard derivation of the virial theorem and obtain, by use of 
the equation of motion (|4.1U|) , 



We may now define the effective frequency ftk as 



Re 



Re 



-fc(^ 3 )fe 



\<h 



(4.11) 



where, with the exception of uJf. 1 everything depends on time, but only according to the slow motion. 

Our aim here is to derive a quasi-particle description for the classical evolution. Quasi-particle pictures have been 
succesfully derived in the context of the Hartree approximation both classically and quantum mechanically (see for 
example refs. [1 H E3 El 13 ) . 
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FIG. 26: The space average of <f>(x,t) as a function of time t for E/V = 100 with L — 6.4 and a — 0.1 over a sample time 
interval. Notice the fast oscillations displayed here which are erased by the time averaging. Thanks to time averaging the slow 
dynamics becomes visible in figs. 121231 



On the staggered lattice there are purely kinematic modification to eq. (|4.11|) . since the fields satisfy discrete 
recursion relations rather than differential equations. In Appendix[B]we derive the virial theorem for the free discrete 
dynamics on the staggered lattice. A comparison of eq. H_B6|) with eq. I|4.11|l the suggests the following lattice relation 
for the effective frequency fife 



l**l a 



1 - (1 - ia 4 ) cos 2 (fi fe a) 



(4.12) 



After averaging over discrete directions, we obtain the following practical rule for calculating the lattice direction- 
independent effective frequency 



(1- ± a 4 )cos 2 (fi fe a)~l-a 2 ^ 



(4.13) 



Eq. H4.13fl suggests also a way to parametrize the equilibrium two-point function (|0fc| 2 ) on the lattice in terms of 

an equilibrium effective frequency fi fc q : replacing |7Tfc| 2 by its equilibrium counterpart (|7Tfe| 2 ) = T and \<fik\ 2 by (\<fik\ 2 ) 
yields 



(10 



a 2 T 



l-(l-i a *) cos2(fi^a) 

Going back to the continuum, we recall that the Hartree approximation for fife follows from the reduction of 



(4.14) 



qq> 



^h—q—q' 



to the sum of all terms proportional to <j>k, that is 



>k - T7 



q 
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FIG. 27: M* s (t) - 1 as a function of the logarithm of the time t for E/V = 0.1, 10, 100 and 1000 with L = 6.4 and a = 0.1 . 
No average is performed over initial packet amplitudes or positions. Time averaging is performed here as by eas. 12. 171 1- 12.1811 . 
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FIG. 28: R = [M*g(t) - l\/4> 2 {t) as a function of the logarithm of the time t for E/V = 0.1, 10, 100 and E/V = 1000 with 
L — 6.4 and a = 0.1. Notice that R varies in a much narrower interval than M^ B (t) depicted in fig. 1271 No average is performed 
over initial packet amplitudes or positions. Time averaging is performed here as by eas. 12.17H - l2.18ll . 



37 



Using ea. H4.11fl this yields the Hartree frequency, 

(0*) 2 = ^+3^(t). (4.15) 

It must be noticed that the same effective frequency follows by using the Whitham approach]^. That is, considering 
a multi-wave configuration as in ea. (|4.1|) . the first nonlinear correction to the frequency in the Whitham approach 
can be shown after some work to coincide with the Hartree formula ea. (|4.15|l . 

Taking into account the definition of time-dependent wavefunction renormalization Z k (t) as defined by ea. l|4.8[l . we 
now have 



Z fe (i)^4=fl-(l-K)cos 2 («^)l =§S£ (4.16) 



a 2 |7T fc | 



Hence 1/Zf.(t) is the renormalization which turns the approximated Hartree frequency into the exact effective fre- 
quency. In the previous section we have provided numerical evidence showing that Zk(t) takes a form quite close to 
the equilibrium one (at some time-dependent effective temperature) for t > to and in particular that it approaches 
unity rather fast as k 2 > cf> 2 (t) . As shown in section IIV El and in more detail in appendix El this holds for very 
late times of order IOHq, close to complete lattice thermalization (necessarily cutoff-dependent), with Zf.{t) almost 
constant in time much earlier than <j) 2 {t) or the power spectra and the effective temperature very close to the final 
equilibrium temperature. But it also holds for much shorter times, of order to only, as soon as the universal cascade 
has set in, with Zk(t) almost constant as a function of fc[(/) 2 (t)] 1 / 2 and some time-dependent effective temperature 
much larger than the final equilibrium temperature. 

We tested the universality of this picture by studying the changes of Z k (or absence thereof) upon changes of the 
lattice spacing, physical size, energy density, initial condensate and microscopic details of initial conditions in our 
numerical evolutions. 

Notice that we always use the exact classical time evolution for the fields and never the Hartree approximation to 
it. However, this exact evolution of the modes with k 2 > 4> 2 (t) is well reproduced with an effective mass as given by 
the Hartree approximation ea. H4.15|l . 

In quantum theory, as is well known the definition of the number operator is not unique (see, for example refs.0,|^]). 
However, we can introduce a classical number of modes based in the correspondence principle. In the classical limit 
the number of quanta is given by the phase space area encircled by the classical phase space trajectory divided by 
2 7r (in units where H — 1). In the Hartree approximation both TTk(t) and <pk{t) oscillate with frequency fl^ (t). The 
phase space area for a slow varying frequency is then given by the ellipse area it \TTk\max |^fc| max where in addition 
for harmonic oscillations |7Tfc| ma:E = (2 |7Tfc | 2 ) x / 2 and |</>fc| m a£c = (2 l^fcl 2 ) 1 ^ 2 - Therefore, the number of modes is given 
by 

1/2 



n fc (t)= |^| 2 (i)|0 fe | 2 (t) . (4.17) 



Also, 



«»(') = !|fw(0 + l||w«) 



where we used the continuum limit of eq. (|4.16(l , namely, 



Z k (t)\n k \ 2 (t)=[nl{t)Y \4> k \ 2 (t). (4.18) 

Although derived within the assumption of harmonic oscillations with a slowly varying frequency, eq. I|4.17|l should 
have a more general validity, since n k provides in any case a measure of the phase space occupied by the trajectory of 
the k mode. The only important condition is that the different modes are weakly coupled. This is true, when the UV 
cascade has fully developed, for all k such that [</> 2 (f)] 1//2 < \k\ < k(t), where [(^(i)] 1 / 2 decreases with time tending 
to a small equilibrium value yj ((f) 2 ) of order y/T '/a ~ a \J E/V. Furthermore, if t is not too large so that k(t) <C A, 
then the occupied modes do not feel the lattice discretization and have the relativistic dispersion relation. Therefore 



n k (t) - l**l 2 (*) > ( 4 - 19 ) 
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where now [n%(t)] 2 = k 2 + 1 + 3 2 (t). 

We plot in fig. EH the number of modes over spherical shells, 47rfc 2 77fc(£), vs. k = \k\ at different times. Fig. |2H1 
should be compared with figs. ITTI and IT51 since nfc(t) and |7Tfc| 2 (i) are related through ea. l|4.f 9|l . One sees that the 
main dependence of rTk(t) on k comes from |7Tfc| 2 (i). 

Notice that eq. 14.19|l can be written in terms of the effective frequency fi&(t) as 



- (A l^l 2 W 

nk{t) = im 



which is the classical equilibrium occupation number with the temperature replaced by |7Tfc| 2 (i). In fact, in thermal 
equilibrium the power spectrum and the temperature are related by eg. 13. 411 . During the UV cascade we are in a 
situation of effe ctive equilibration for times t later than to ~ 500 as shown by fig. |5]and the results of section flV Fl 
However, | tt/c | 2 (i) depends both on time and on the wavenumber as depicted in figs. ITU HTl and l2~3l In addition, the 
k- modes with k 2 < (j) 2 (t) only thermalize for times of the order 10 6 as we see in fig. [23 This makes it awkward to try 
and interpret |7Tfc| 2 (i) as a k— dependent effective temperature since an equilibrium or quasi-equilibrium state should 
depend only on few macroscopic parameters (with one playing the role of temperature) varying slowly in time. Notice 
however that |7Tfc| 2 (t) decreases both with time and with the wavenumber k for k 2 > 4> 2 (t) [see figs. If 41 If 61 and!23 j 
and does that in a very smooth and regular way. This suggests indeed the existence of few slowly time-dependent 
parameters, not dependent on th e details of the initial conditions, which govern the evolution of the cascade for 
k 2 > 4> 2 (t). On the contrary |7Tfe| 2 (t) increases as well as decreases with k and t for the modes k 2 < 4> 2 (t) in a much 
more erratic way strongy dependent on the detail of the initial conditions [see figs. If 41 1221 and I23 j . 



H. Effective Mass squared 

The time averaging defined by ea. (|2.f 7|) is intended to eliminate the microscopic oscillations of the field. We 
illustrate such microscopic behaviour in fig. 1261 showing the field <fi(x,t) averaged over the space. The frequency 
of the time oscillation can be therefore considered as the effective mass M g{t) of the field. Moreover, by Fourier 
transforming the field 4>{x, t) we obtain frequencies fife numerically. This procedure is more costly than eq. (|4.f 3|) 
but we find an excellent agreement between both. Notice that the use of lattice expressions like eq. (|4.f 3f) . valid in 
principle to all orders in a, is very efficient and convenient to compare results on lattices with different lattice spacings. 

It must be noticed that the fast oscillations displayed in fig. [201 are erased by the time averaging. Thanks to such 
averaging the slow dynamics becomes visible through figs. 121231 

We plot in fig. |23 [M 2 ff (i) - I] vs. the logarithm of time for different values of p = E/V. Notice that M 2 s (t) 
monotonically decreases with time while the UV cascade develops. We see that M 2 S (t) > I and that it decreases with 
p. Indeed, M 2 s (t) — > 1 (its linearized value) for p — ► 0. We find that M 2 s (t) — I is approximately proportional to 

We depict in fig. [2H1 the ratio 

R= M ^-\ (4.20) 

as a function of \ogt. We see that R is approximately time independent and that it decreases with E/V. It should 
be noticed that R always stays below the value Rnartree = 3 predicted by the Hartree approximation [|| and well 
above the value RLargeN = I corresponding to the large number of components limit. 

In summary, we find that M 2 s (t) can be written as 

M 2 s {t) = l + R{p)^{t) , (4.21) 

where 1.5 < R{p) < 3 is approximately time independent and decreases with p. 

We analyze the cnoidal solution in appendix C and compute its effective mass. The parameter R for the cnoidal 
solution turns to be between | and | and is depicted in fig. El 

As discussed in sec. II V Gl the modes with k 2 3> 4> 2 {t) are governed by the Hartree mass I + 3 <j) 2 (t) while the k = 
mode oscillates according to M e g(t). The Hartree mass and M 2 s (t) are substantially different as remarked above. 

Moreover, M c g (t) governs the long range behaviour of the correlators is real space. The two-point correlator tends 
to zero for long distances as, 

„ p-Afeff(t) X 



C{x) = Co 



x 
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FIG. 29: The number of modes over spherical shells, 4-7T k 2 nfc(t), vs. fc = |fe| at the different times indicated. The initial 
conditions are infrared random plane waves. The UV cascade is clearly seen as k 2 rifc(t) is depleted for low k while it grows for 
large k as time grows. 

where Co is a constant. 

In summary, modes with k 2 below and above (j> 2 (t) are governed by different masses. The first by M 2 s (t) ~ 
1 + 1.5 <j> 2 (t) and the later by the Hartree mass 1 + 3 4> 2 (t). Clearly, the long-range behaviour of the correlators is 
governed by the low momentum mass M e g(t). 

I. The slow dynamics of the infrared modes and effective thermalization 

The picture of the averaged dynamics for late times in the lattice model is as follows: 

• The modes with k 2 > 4> 2 {t) effectively thermalize with a time dependent temperature while they obey the 
Hartree approximation. That is, their dynamics is weakly nonlinear. 

• The modes with < k 2 < 4> 2 (t) keep interacting quite strongly between them. They thermalize much later than 
the rest of the modes 

• Both set of modes keep interacting between each other 

There are definitely two time scales for thermalization. The shorter one to ~ 500 describes the effective thermaliza- 
tion of the Hartree modes while time scale for the effective thermalization of the infrared modes is much longer than 
to and of the order ~ 10 6 . 

The modes with < k 2 < (j> 2 (t) manifest in the IR behaviour of the correlator as a long-lived inhomogeneous 
condensate. This condensate keeps interacting with the higher k modes which behave as a thermal bath in contact 
with it. Notice that the mode distribution becomes time independent on the lattice for late times. Hence, the field 
correlator is a static function described by the Fourier transform of 
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where T is the asymptotic time limit of |7ffc| 2 (i). 

While Mg ff decreases with a time scale ~ to [see fig. EZ], %k varies with an even longer time scale ~ 10 6 . It can be 
interpreted as a wavefunction renormalization which becomes unit for k 2 ^> (f) 2 (t) . 

It should be stressd that the inhomogeneous condensate is not a classical solution of the <f) A equations 1)2. 3)1 . It is a 
classical statistical configuration that can only be seen through the correlators like (f><f)(x,t). Moreover, the average 
value of (p just vanishes for late times as shown in fig. [H] 

Notice that the local observables as 7r 2 (£), (V0) 2 (i), 4> 2 {t) and <j>+(t) are dominated by the modes with k 2 > 4> 2 (t). 
This explains that the thermal equilibrium curve approximately agrees (to 90%) with time averages in fig. [H] although 
the infrared modes are not yet completely thermalized. 

Extensive numerical calculations varying the initial conditions showed us that the thermalization dynamics including 
the UV cascade posseses an universal character. That is, the UV cascade exhibit the same features for different kinds of 
initial conditions (plane waves, localized packets, etc.) with initial power in the infrared modes (k < 30 -k / L <C t/[2 a]) 
and energy density not too small (p > 1). 

In the continuum theory the ultraviolet cascade continues forever and therefore, for finite energy densities, the 
system will reach thermal equilibrium for infinite times at zero temperature. The inhomogeneous condensate thus 
disappears for extremely long times in the continuum theory. 

As the inhomogeneous condensate is an infrared phenomena, it should be present in quantum theory for large 
occupation numbers and weak coupling. In other words, the low k behaviour of Z% should not change in quantum 
theory when such conditions are fulfilled. 



V. DISCUSSIONS AND CONCLUSIONS 



The present work shows that a UV cascade enjoying universal properties is the basic mechanism of thermalization 
both in 3 + 1 and 1 + 1 dimensions 6] . Although we are not in a position to derive the properties of this cascade 
from the microscopic equation of motion 1)2.3)) . it is easy to see from ea. <)2.3)) that energy should flow towards higher 
wavenumbers. Assume we start from a plane wave initial condition 

0(0, x) = A cos(fc • x ) , <j)(0, x) — A sin(fc • x) , 

where A is a constant. The nonlinear term (j) 3 in Eq. (|2.3(l will immediately generate higher harmonics: cos(3 k ■ 
x), cos(9 k ■ x), cos(5 k ■ x) etc. That is, energy is moving to higher wavenumbers. 
Moreover, if we consider a superposition of plane waves as initial condition, 

n n 

4>(0,x) = ^ A a cos(k a ■ x) , 0(0, x) = A a sin(fc a ■ x) , 

a— 1 a— 1 

the interaction term <fi 3 generates sums and differences: k a ± kf, ± k c etc, for a, b, c = 1, . . . , n. This implies energy 
flowing both for increasing and for decreasing wavenumbers. 

Notice that this mode mixing takes place very fast, at the microscopic time scale. That is, a unit time scale in 
dimensionless variables. The UV cascade we observe happens in a much longer time scale. This means that these fast 
microscopic processes combine in a nontrivial way resulting on a UV cascade with a slow time scale. 

It is important to estimate in which extent the present results can be applied in quantum theory. A necessary 
condition for the validity of the classical approximation is that the occupation numbers must be large. The relevant 
occupation numbers decrease during the UV cascade due to the modes flow towards unoccupied high wavenumber 
slots. Hence, if the classical approximation is valid for late times, it will be valid earlier (as it is the case 1 + 1 
dimensions J5]). A simple criterion for the validity of the classical approximation in QFT at thermal equilibrium is 
that 

> Wp(fcp) (5.1) 



where Lo p (k p ) — Jk^ + M 2 ,? are the dimension-full frequencies. 
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In dimensionless variables eg . (|5 . 1|) takes the form, 

j » ^fc 2 + M c 2 ff > M eff (5.2) 

Following ca. (|4.21(l we obtain as an estimate for M 2 ff from our numerical results 

M c 2 ff ~ 1 + A VP (5.3) 

where p = y and A ~ 0.2 - 0.5. 

For p > 1 we thus find that the classical approximation applies for 

A < (2a) 3 pi (5.4) 

where we used that T = (2a) 3 p for small a [ea.i|3.9|)]. 

For p < 1 we instead find that the classical approximation applies for 

A<(2a) 3 p (5.5) 

Both eq. (|5.4|) and 15.511 leads to small values of the coupling A since the spacing a must be itself small to avoid lattice 
effects. However, in inflationary theories very small values of A ~ 10~ 12 are customary (in order to agree with the 
smallness of the CMB anisotropy) which leaves room for the use of the classical approximation. 

As stated above the validity of the classical approximation decreases with time. Let us estimate the time where it 
ceases to be valid. 

The total energy density can be estimated in terms of the average occupation number (where the average is here 
over the k modes) as follows 



3 / 2 

p ~ k rik \Jk + M 2 ff . 

The classical approximation holds if the occupation numbers are large for the relevant modes. Therefore, a necessary 
condition is 

nk 

-T > 1 i 
A 

since the coupling A has been absorbed in the field redefinition Eq. I|2.2(l . 
In the ultrarelativistic regime k 3> M f{ we have 

k 4 < j and therefore, k < (J^J 4 , (5.6) 

while in the non-relativistic domain k -C M c q we have, 



and hence using ea. l|5.3|l 



Since A <C 1, 



, r t3 P 

M eS k < j- 



—3 pi —pi 

k -C — and therefore, k <C — r ■ (5-7) 
A As 



1 1 , { P\i pi 

-r > -r , thus ( - < -r , 
A3 \i V A/ As 

and the ultrarelativistic condition eq. (|5.6(l is more stringent than the non-relativistic bound eq. (|5.7() . We can therefore 
always use the condition eq. (|5.6|l for the validity of the classical approximation. 

Using now ea. l|4.7|) for k(t) yields in both regimes that the classical approximation is valid for times t 

3 

t <C t max ~ —3 — - . (5-8) 
A* kR 
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The conditions eq. (|5.5() and l|5.6|) are compatible since k < ^ . Therefore, eq. I|5.6[) implies 

7T- < It i * and tncn ( 2a ) 3 P^ > 71-3 or t > 
2a V A/ A 



-V 

2a/ 



which is more stringent than ea. H5.4|) . 

In summary, the classical approximation is valid for times earlier than t max [given by ea. l|5.8|l ] provided 



p / IT \ 4 



That is, high density (large occupation numbers) and/or small coupling. This condition constraints the initial 
conditions which fix p. Notice that the coupling A does have an intrinsic meaning in the quantum theory. 

A condition for the validity of the classical dynamics in QFT is derived for high-density in ref. in the context 
of the 2PI-1/N approach. For further studies about the validity of the classical approximation in different contexts 
see [liill^. 

Let us now comment about the character of the universal stage. The dimensionality of the space plays a crucial role 
in this phenomena. In one space dimension effective thermalization takes place much faster in analogous conditions Q. 
Here, in 3 + 1 dimensions a second and even longer time scale appears characterizing the thermalization of infrared 

modes with k < \J~cjp(t). The infrared modes keep interacting for far longer times than the higher k modes. In 3 + f 

dimensions the infrared modes thermalize the last while in 1 + 1 dimensions they thermalize the first 0. 

The UV cascade is clearly less efficient in three space dimensions to fill the higher fc-modes pumping energy from 
the low fc-modes. Clearly, in one space dimension the phase space is dramatically small making the UV cascade very 
efficient. 

The same phase space effect (k 2 ) in three space dimensions makes the classical statistical mechanics of the (/> 4 theory 
divergent due to the ultraviolet catastrophe. 

All this suggest that thermalization is reached in quantum theory too as recent works (including memory effects) 
indicate |l4l Il5l IT^ . Moreover, for large initial occupation numbers and small coupling the classical regime should 
correctly describe the evolution of the theory till a time t ma x [see ea. (j5.8(l ]. 

A remarkable feature of the thermalization mechanism is that even starting from a completely classical regime (large 
occupation numbers at low wavenumbers), the UV cascade depletes these modes and fills the high ones. Therefore, 
at some time ~ t max (that could be very long) quantum physics unavoidable shows up. This is the out of equilibrium 
counterpart to the fact that classical statistical mechanics is ill defined due to the ultraviolet catastrophe which can 
only be cured by the quantum treatment. 

Thermalization implies to forget everything about the initial conditions except (obviously) the conserved quantities 
like energy, momentum, angular momentum. The thermalization is therefore expected to substantially change for 
integrable theories where the number of conserved quantities equals the number of degrees of freedom. 

In this paper we solve the exact microscopic dynamics and then we averaged on time intervals and space in order 
to derive the slow dynamics we are interested in. This is perfectly correct but a lot of information is first obtained 
and then dumped in the averaging process. Alternatively a transport approach could be followed deriving equations 
for the observables in macroscopic time scales and then solve them. We mean to derive transport equations of the 
Boltzmann or Fokker-Planck type. 
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APPENDIX A: AVERAGE OVER DISCRETE DIRECTIONS 



The space discretization and the finite size spoil the rotational invariance of continuum infinite-volume Hamiltonian. 
This implies the same invariance for the classical field equation ad for the equilibrium averages (see section llllll . 
Rotational invariance should be recovered at short distances when a — ► and at long distances when L — > 00. Thus 
any lattice observable on the wavenumbers cube (2n/L) C[/v that corresponds to a rotational invariant observable 
of the continuum infinite-volume theory, should depend only on k 2 when a — > and kL 3> I or when L — > 00 and 
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ka <C 1 (a similar argument would apply for the x — space cube 2a Cat)- To verify the onset of rotational invariance 
in non-perturbative or numerical lattice calculations is usually quite costly and often not even necessary. In fact, by 
turning the argument around and assuming that rotational invariance will be recovered in the appropriate limits, one 
can greatly reduce fluctuations by performing an average over all directions on discrete observables. This can be done 
as follows. 

Given any array f n over Cat, we consider its average over discrete directions f n as 

fn = ^~ J2 0(n<\n\<n+l)f n , n = 0, 1,2, . . ., V3iV/2 . (Al) 

where 0{. . .) = 1 (0) if its arguments is true (false) and S n is the number of points of Cat at a distance d from the 
origin such that n < d < n + 1 , that is 

S n = VJ 9(n< \n\ <n+l) 

Up to purely geometrical fluctuations, S n grows like 4 tt r? for n up to N/2, while for N/2 < n < \/3 N/2, in the 
corners of the cube, its shrinks to zero. To fix a specific value of the radius of this spherical shells, we choose the 
average distance of all lattice points within the shell 

r„ = ^- 0(« < M < n + 1) \n\ 

When plotted against r n , S n still has purely geometrical fluctuations w.r.t. the continuum expression, most noticeably 
at n = 0, where Sq — 1 while Airr^ — 0. Of course, averaging over several consecutive shells would reduce these 
geometric fluctuations. Such an average is in fact automatic in the limit L — > oo of continuous wavenumbers, since 
in any physical realizable measure there exists a finite resolution Ak independent of L. On the other hand the limit 
L — ► oo is not feasible in numerical calculations and one has to find the correct balance between size and control of 
fluctuations. 

Suppose now that the array /„. if the finite size version of the continuous function f(k) over the first Brillouin zone 
at infinite size (the argument would apply also if f n were the lattice version of a continuous fix) on the cubic volume 
V). Then /„ provides a finite size version of the angle average of f(k) of f(k), where k — \k\. In other words, the 
discrete plot of /„ vs. k n = (2n/L) r n provides an approximation of the continuous plot of f(k) vs. k, if L — > oo at 
fixed a and f n has indeed a continuous limit. If this limit is rotational invariant all physical information is stored in 
/ as a function of k only; thus the average over all solid angle in the continuum has no effect at all, f — f, while the 
average over discrete shells S n greatly reduce the statistical fluctuations. 



APPENDIX B: LINEARIZED DISCRETE DYNAMICS 



It is instructive to study first the discrete field equation (|2.12|) at the linearized level. Thus we consider the case of 
a uniform field <fi(x,t) — 4>(t), namely the recursion 



2 0(0 
1 + | a 2 [l + <f{t)} 



4>{t + a) + ^{t-a) = - ; l J) \ (Bl) 



Given 0(0) and </>(a), this recursion determines cj>(t) for the discrete times t = na, n = 2,3, . . ., yielding the discrete 
version of the well-known cnoidal uniform solution of the continuum equation (|2.3I) . 

We now linearize the generic 4>{x,t) around cf>(t) by setting 4>(x,t) — <fr(t) + r)(x,t), where r){x,t) has vanishing 
space integral; this yields the linear recursion 

r)(x, t + a)+ r)(x, t - a) = - B((f>) r)(x + acr, t) (B2) 

<T 

where B(cj)) is given by eci. (|3.2(j|) In terms of Fourier modes we have 



r,{x,t) = V- 1 ' 2 



ik-x 
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and eq. I|B2(1 entails, for k ^ 0, 

4>k(t + a) + 4>k(t - a) = 2 4> k (t) coswfc(t) 
where the instantaneous frequency uj k (t) is fixed by the dispersion rule 

D 

cos [ujk(t) a] = B[4>(t)] Y[ cosfcja 



(B3) 



(B4) 



as function of the wavenumber vector k ^ 0. Eq. (|B3|) is our lattice version of Lame equation governing the 
linearization of the continuum field equation over the cnoidal. In fact, in the continuum limit a — > 0, for any 
fixed fc we have 



cos [w h (t) a] = 1 - ±a 2 [fc 2 + 1 + 3 <f> 2 {t)} + C(a 4 ) 



so that eq. I|B2|) becomes indeed 



dt 2 



+ fc 2 + l + 3 2 (t) 4> k (t) = o 



(B5) 



with (f>(t) solving now <f> + (j> + <fi 3 = 0, that is the cnoidal solution. 

In the trivial case <f> = 0, which corresponds to the free field case, the general solution of eq. (|B3(I is straightforward 
and can be written 

4> k (t) = A k e-^ kt + A*_ k e ZUJht 

in terms of the constant amplitudes A k fixed by the initial conditions. In this case u> k solves eq. <|B4(1 with B(0) = 
(1 + a 2 /2) _1 and is constant in time. Then recalling the general relation eq. I|3.22|l we obtain 



7Tfe(*) 



.sin (wfcd) 



(A k e-' 1 ^ - A*_ k e^*) =p y &(t) cos (w k a) 



Finally, since time averages over several mode oscillations kill the interference between the positive frequency and 
the negative frequency terms in the power spectra, we find in this <p = case the following lattice version of the well 
known virial theorem for harmonic oscillations 



W = 



l[l-(l-ia 4 ) cos 2 Ka)j |0 fc p = J_h_:L4^ f[ 



cos 2 fco-a 



(B6) 



2- j= l 

where the dispersion relation eq. (|B4(1 was used with (f> = 0. In the limit a — > eq. (|B6|) reduces to 

= (fc 2 + 1) j^(t) 

as expected. 

When ^ we may still write the general solution of the recursion (|B2(I as 

Mt)=A k z k {t)+A*_ k z* k (t) 
where the amplitudes A k are constant in time while z k (t) = z^ k (t) solve eq. <|B2|) with initial conditions 

z k (0) = 1 , z k (a) = e -~*C°)« . 



Eq. I|B5|I admits close form solutions when </> is given by the cnoidal solution ea. (|Dl(l . In particular, the forbidden 
band for fc 2 > corresponds to 



\ $ + 3 < fc 2 < ^$ + 2^ + 4+1 



(B7) 



(see for example [Til). 
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APPENDIX C: THE HARTREE APPROXIMATION FROM THE LATE TIME EXACT BEHAVIOUR 

We learn how to obtain the Hartree frequency on the lattice: we consider the discrete dynamics linearized as in 
Appendix iBl and replace the background </>g(t) with <fi 2 (t) in the dispersion relation eq. (|B4|I . obtaining 



cos(af4 H ) = B[£0 2 (i)] JJcosfya 



(CI) 



where -B(</>) is given by Ea. (|3.2t)|) and where for future use we introduced the parameter £ to control the coupling 
with the background: if £ = 1 we have the Hartree frequency, if £ = we have the free massive frequency ea. H4.91 . 
Intermediate values < £ < 1 will actually be ruled out below by fitting the numerical data. 

f^ H is a function of time through the background (f> 2 (t). On the lattice (j> 2 (t) will tend to a nonzero limit as 
t — > oo. According to the Hartree dominance of the equilibrium (\<fik\ 2 ), from eq. i|4.14[l we than see that f2fc ~ 57^, 
provided the full effective frequency flk tend to £!^, q as t — > oo. It is therefore very important to compare £!^, H 
with the full effective frequency fife. Recall that the Hartree approximation is a good approximation if the modes 
are weakly interacting, since it neglects direct wave scatterings and allows energy transfer only through the uniform 
self-consistent background cf> 2 (t). 

To perform the comparison we numerically average over discrete directions and study the ratio 



AO 



a 2 |fr fc | 2 



I - (I - ±a 4 ) cos 2 (a fif) 



(C2) 



In fig. 123 we plot Zj^ , Zj^' 5 ' , zj^ and Z^' 5 ^ for several late times in the long evolution when ir thermalized (see 
previous section). The Hartree value £ = 1 clearly stands out as the most accurate in the time dependence: the curves 
at all different times almost perfectly collapse in a single curve, except that for very small wavenumbers. All curves 
collapse for large k, regardless of the value of £, simply because when k 2 ^> <fi 2 (t) any effect on the dispersion relation 
due to the background is negligible. 




1 = 2242.82 
3187.33 
5351.04 
11464.7 
30376.4 
90753.7 
285479 
915494 



A = 24.543 



L = 12.8 
a = 0.064 
E/V = 569.5 



10 



15 



20 



25 



FIG. 30: The inverse of the ratio Z\ [defined by eq. HC2H vs. the radial wavenumber for the values of time and of £ indicated. 
All parameters are as in fig. 1201 
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APPENDIX D: THE CNOIDAL SOLUTION 



We discuss here the cnoidal solution of the cf) 4 equation (|2.3|l . that is the space- independent solution of 4>+<j)+ (j> 3 = 
explicitly given by 



4>(t) = 4> cn^t \Jl + 4>1, kj with k 



V2(l + ^) ' 

where cn(x, k) stands for the Jacobi cosinus of module k. This is a (doubly) periodic function satisfying <f>{t) 
—4>(t + P), where P = 2 K(k) /yl + </>q is the half-period and K(k) stands for the complete elliptic integral. 
The effective mass associated to this solution is therefore the basic frequency 



(Dl) 



A I, 



off 



2_7T _ 7T yj\ + (j)l 

2P 



2 K{k) 

The average value of the squared field over a period can be computed with the result (see eq. (3.40) in ref.||). 



1 

2P Jo 



2 P 



dttf(t) = {l + cj>l) 



2 E(k) 
K(k) 



- 1 



(D2) 



(D3) 



We plot in fig. EUthc ratio R(<f) ) = (Af c 2 ff - l)/((f> 2 ) using Eqs. (JD2|) and |jD3l) [compare with ca. (|4~2TJ|) ]. Wc plot it 
as a function of the energy density p — \ 0§ (l + \ </>q) . 

We can compute analytically the function R(4>q) for small and large arguments with the result 



i>o^o 3 



In summary, R(4>q) grows slowly and monotonically from the value | to the value ^ when (j>o varies from zero to 
infinity. 

The average value of 4> 2 can be computed analogously with the result (see eq. (3.60) in ref. 9]), 



2 > = 3(1 + ^0) 



We find from Eqs. (|D3|) and l|D4|) using energy conservation 



2 - 



1 



2 E(k) 
K{k) 

E{k) 



3 K{k) 



(1 + Po) 



(D4) 



(D5) 



Inserting Eqs. (|D3(I . I|D4|I and l|D5(l for this homogeneous solution in the virial theorem Ea. H3.6(l shows that this 
theorem is identically satisfied. 
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